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A  nurr,erical  model  is  employed  to  describe  the  water  motion  in  a 
coastal  region  associated  with  the  passage  of  a  hurricane  or  severe 
storm.  The  model  is  two-dimensional,  employs  the  vertically  integrated 
or  "tidal"  equations  of  motion,  and  is  used  to  describe  the  specific 
case  of  Hurricane  Camille  of  August  1969.  Two  models  are  employed,  a 
large  grid  with  16  nautical  mile  grid  elements  and  a  small  grid  model 
with  6  nautical  mile  grid  elements.  The  results  of  the  two  models  were 
not  significantly  different.  The  maximum  storm  surge  or  maximum 
water  level  above  mean  sea  level  was  calculated  and  found  to  agree  to 
within  twenty  per  cent  with  maximum  surge  heights  determined  from 
high  water  marks.  Calculated  surge  height  was  found  to  be  insensitive 
to  bottom  friction  coefficients  varying  from  .005  to  .02.   Inclusion 


XV 


of  the  nonlinear  convective  terins  affected  the  calculated  surge 
height  maximum  for  Hurricane  Camille  by  less  than  two  per  cent  and 
affected  the  surge  in  any  grid  element  adjacent  to  the  coast  by  less 
than  five  per  cent.  Tv/o-dimensional  storm  surge  plots  at  different 
times  are  presented. 

Hurricane  generated  currents  were  calculated  and  compared  to  data 
taken  off  the  Florida  coast.  It  is  concluded  that  more  actual  current 
data  are  necessary  before  hurricane  generated  currents  can  be  cal- 
culated with  confidence.  The  current  calculations  were  found  to  be 
sensitive  to  bottom  friction  and  subject  to  "wind-up". 

In  addition  to  numerical  calculations,  basic  analytical  cases 
are   covered,  including  the  response  of  a  shelf  of  uniform  depth  to 
a  triangular  wind  stress  distribution  moving  with  constant  velocity. 
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CHAPTER  1 

INTRODUCTION 

The  problem  of  predicting  the  response  of  coastal  waters  to 
extreme  wind  systems  is  becoming  increasingly  important  as  the  coastal 
regions  of  this  country  are  developed  for  residential,  recreational 
and  industrial  uses.  Disasters  such  as  the  ones  occurring  in  recent 
years  at  Galveston  as  a  result  of  Hurricane  Carl  a  and  along  the 
Mississippi  coast  as  a  result  of  Hurricane  Camille  demonstrate  that 
action  should  be  taken  to  minimize  the  tragic  consequences  of  such 
events.  Prediction  of  the  vulnerability  of  particular  areas  to  flood- 
ing in  similar  storms  would  be  of  considerable  interest  to  local 
residents  and  governments  as  well  as  insurance  underwriters.  Moreover, 
real  time  forecasting  of  storm  surge  elevation,  with  confidence,  for 
an  approaching  storm  would  be  of  great  value  in  planning  and  imple- 
menting evacuation  and  other  measures  to  minimize  damage.  This  is 
clearly  illustrated  by  Hurricane  Agnes  (1972).  Agnes  was  a  hurricane 
of  large  areal  extent  but  relatively  low  intensity.  Little  attention 
was  paid  to  surge  forecasting,  and  as  a  result  the  residents  of  some 
areas  on  the  west  coast  of  Florida  were  surprised  by  a  storm  tide  3  to 
8  feet  above  normal.  Damage  in  the  Tampa-St.  Petersburg  area  was 
estimated  at  about  $20  million.  Undoubtedly  some  of  this  could  have 
been  prevented  if  adequate  real  time  surge  forecasts  had  been  avail- 
able. 


The  design  of  permanent  facilities  at  or  near  the  shoreline 
must  include  consideration  of  the  maximum  and  minimum  tides  that  may 
occur  as  well  as  the  wave  action.  In  the  operation  of  nuclear  power 
plants,  for  example,  one  consideration  is  the  availability  of 
sufficient  cooling  water  to  ensure  an  orderly  shut  down  prior  to  a  low 
tide  occurence  below  the  level  necessary  for  plant  operation.  Thus, 
a  proper  study  of  the  situation  would  establish  a  maximum  elevation  at 
which  the  inlet  facilities  could  be  placed. 

Various  approaches  are  used  in  solving  this  type  of  problem. 
The  purely  analytical  approaches  are  limited  to  idealized  wind 
stresses,  pressure  distributions,  and  topography.  Lamb  [1]  considers 
the  case  of  "free"  long  waves  in  an  infinite  canal  of  uniform  depth. 
In  addition  to  gravity,  small  disturbing  forces  are  considered  to  be 
acting  on  the  fluid,  which  vary  only  by  a  small  fraction  of  their  value 
within  distances  comparable  to  the  depth.  This  leads  to  an  in- 
homogeneous  wave  equation  similar  to  that  of  Section  3-e  of  the  text. 
Lamb  then  finds  a  solution  for  a  pressure  distribution  that  travels 
with  unchanged  form  at  constant  velocity.  The  presence  of  boundary  or 
terminal  conditions  greatly  complicates  the  solution.  The  case  of  a 
finite  shelf  is  considered  in  the  text. 

Reid  [2]  presents  the  approximate  response  of  a  sloping  shelf 
due  to  a  triangular  shear  stress  distribution  traveling  directly  on- 
shore. The  results,  obtained  at  the  shoreline,  are  based  upon  the 
linear  one-dimensional  wave  equation.  The  solutions  are  obtained 
using  a  graphical  technique  and  the  method  of  characteristics.  He 
considers  many  different  cases  of  fetch  lengths  and  storm  speed  and 
summarizes  his  results  in  graphical  and  tabular  form. 


Further  analytical  work  in  the  area  of  storm  surges  has  been 
and  will  be  done,  yet  the  complexities  of  the  actual  forcing  functions 
and  coastal  topographies  lead  one  to  realize  that  another  approach  is 
necessary,  if  realistic  solutions  to  actual  problems  are  to  be 
obtained.  Hansen  [3,4]  is  one  of  the  pioneers  of  storm  surge  cal- 
culations. He  calculated  surges  arising  from  storms  occurring  in  the 
North  Sea,  using  the  linearized  equations  of  motion.  That  is,  he 
ignored  surge  height  with  respect  to  depth  and  ignored  the  convective 
terms.  He  laid  the  foundation,  however,  for  the  present  work  in 
numerical  surge  calculations.  Many  contemporary  investigators  still 
use,  in  some  form,  the  staggered  integration  scheme  he  devised. 

Platzman  [5]  analyzes  the  case  of  a  pressure  distribution  or 
squall  line  moving  over  Lake  Michigan  causing  a  surge  of  more  than 
six  feet  on  the  Chicago  lake  front.  The  volume  transport  or  ver- 
tically integrated  equations  of  motion  were  used  in  central  difference 
form  in  conjunction  with  a  two-dimensional  bathymetry  and  a  moving 
pressure  disturbance.  Platzman's  calculated  surge  is  about  one- 
half  of  the  observed  surge  and  he  concludes  that  this  is  probably  due 
to  the  omission  of  wind  stress. 

Miyazaki  [6]  treats  the  Gulf  of  Mexico  as  a  closed  basin, 
assuming  the  surge  height  to  be  zero  at  the  mouth  between  Yucatan, 
Cuba  and  Florida.  This  procedure  eliminates  the  difficult  problem 
of  specifying  lateral  boundary  conditions  in  open  water.  The  grid 
he  uses  is  variable  and  finest  at  points  of  interest  in  the  coastal 
regions.  He  includes  bottom  stress  in  a  form  given  by  Reid  [7]  which 
is  a  function  of  both  flux,  q,  and  of  surface  shear  stress,  t  . 

Marinos  and  Woodward  [8]  calculate  storm  surges  using  the 


assumption  of  bathystrophic  flow  (currents  parallel  to  the  bottom 
contours).  They  arrive  at  a  "quasi -one-dimensional"  scheme  which 
considers  only  one-dimensional  depth  profiles  but  includes  the 
Coriolis  effect.  The  bathystrophic  procedure  was  originally  developed 
by  Freeman,  Baer  and  Jung  [9]. 

Reid  and  Bodine  [10]  present  a  detailed  numerical  model  for 
the  prediction  of  surges  in  Galveston  Bay,  Texas.  The  bay  is  treated 
as  a  nearly  closed  system  which,  again,  eliminates  the  problem  of 
open  water,  lateral  boundary  conditions.  Measured  tidal  data  at 
the  entrances  to  the  bay  are  used  as  an  input  to  the  model  and  flow 
through  these  areas  is  allowed.  The  model  includes  bottom  friction 
and  allows  for  over-topping  and  flooding  of  adjacent  low-lying  areas. 

Harris  and  Jelesnianski  [11]  and  Jelesnianski  [12,13,14,15] 
have  presented  a  series  of  papers  concerning  storm  surge  calculations. 
The  numerical  technique  follows  that  of  Platzman.  In  the  paper 
"Numerical  Computations  of  Storm  Surges  with  Bottom  Stress,"  Jel- 
esnianski [14]  includes  bottom  stress  in  his  calculation,  using  an 
eddy  viscosity  and  a  slip  coefficient.  Jelesnianski  notes  that  for 
a  storm  not  moving  too  slowly,  the  bottom  friction  term  is  not  im- 
portant and  may  be  omitted,  but  for  slow  storms  or  storms  moving 
parallel  to  the  coast  bottom  friction  must  be  included.  In  the  paper 
"Bottom  Stress  Time-History  in  Linearized  Equations  of  Motion  for 
Storm  Surges,"  Jelesnianski  [15]  abandons  the  slip  coefficient  and 
maintains  the  eddy  viscosity  with  the  assumption  that  it  is  constant 
in  both  the  vertical  and  horizontal.  Good  agreement  is  obtained 
between  surge  height  data  and  the  calculations.  He  notes,  however,  that 
"currents  are  only  of  casual  interest  compared  to  the  slope  or 


height  variations."  In  general  this  is  true,  yet,  it  is  becoming 
apparent  that  the  currents  are  "^ery   important  in  terms  of  beach 
erosion  and  forces  on  various  ocean  structures.  The  development  of  a 
more  generally  applicable  dissipation  mechanism  would  decrease  the 
sensitivity  of  the  current  calculations  thus  decreasing  the  necessity 
for  "calibrating"  a  model  which  may  mask  errors  and  deficiencies  in 
the  model  by  appropriately  selecting  bottom  friction  coefficients. 

The  volume  transport  equations  are  formulated  in  terms  of 
flux  components  and  do  not  describe  the  velocity  distribution  near  the 
bottom,  which  is  related  to  the  bottom  friction.  Jelesnianski  [15] 
deals  with  the  problem  of  bottom  friction  by  providing  bottom  shear 
stress  as  a  convolution  integral  of  the  surface  shear  stress  function. 
In  this  technique  he  linearizes  the  equation  of  motion  by  making  the 
assumption  that  the  surge  height  is  wery   small  compared  to  the  depth, 
n  <  <  h.  In  doing  so  he  precludes  flooding  in  low-lying  areas  and 
calculations  for  s/ery   shallow  waters.  In  this  report  it  is  antici- 
pated that  the  modeling  of  flooding  in  low-lying  areas  and  surge 
calculations  in  very  shallow  water  will  both  be  desirable,  thus  a 
quadratic  bottom  friction  function  which  allows  the  use  of  the  non- 
linear equations  is  used.  This  is  equivalent  to  making  the  assumption 
that  the  vertical  velocity  profile  is  maintained  in  a  near  steady 
state  condition.  Future  development  of  the  model  will  rely  on  measure- 
ments to  provide  insight  into  the  processes  involved. 

Surprisingly  little  work  has  been  done  in  measuring  hurricane 
phenomena.  The  storm  surge  data  available  are  primarily  from  high 
water  marks.  Tidal  records  are  predominantly  from  enclosed  bays  and 
rivers  making  comparison  to  the  numerical  model  difficult.  The  only 


known  current  data  for  Hurricane  Camille  were  taken  by  Murray  [16,17] 
on  the  Florida  Panhandle.  These  currents,  measured  in  the  surf  zone, 
were  due  to  wave-induced  mass  transport  as  well  as  wind  shear  stress  and 
water  slope.  Clearly  more  data  of  this  type  are  needed. 

The  difficulty  of  reaching  offshore  areas  in  a  hurricane  makes 
manned  observation  of  these  phenomena  nearly  impossible,  so  alternate 
measures  must  be  found.  The  appropriate  instrumentation  must  be 
either  in  place  continuously  waiting  for  a  storm  or  it  must  be  placed 
in  the  path  of  the  storm  prior  to  the  development  of  adverse  sea 
conditions  in  that  area.   In  the  first  case,  the  instrumentation  sites 
must  be  periodically  serviced  with  the  knowledge  that  the  wait  for  a 
hurricane  may  be  a  long  one.  In  both  cases,  there  is  the  problem  of 
instrument  recovery.  An  inexpensive  underwater  gage  and  integrating 
current  meter  have  been  developed  and  are  described  in  the  report. 


CHAPTER  2 
STORM  SURGE  CALCULATIONS  -  THEORY 

a.  Volume  Transport  Equations 

The  equations  of  motion  and  continuity  as  used  in  the  develop- 
ment of  the  storm  surge  model  are  developed  qualitatively  in  this 
section  with  the  details  of  the  integration  over  depth  presented  in 
Appendix  A.  The  coordinate  system  and  nomenclature  are  presented  in 
Figure  2-1.  The  origin  of  the  coordinate  system  is  arbitrarily  chosen 
at  the  shelfbreak  and  is  approximately  centered  in  the  region  of 
interest.  Tfie  x  axis  is  directed  toward  shore,  y  is  oriented 
approximately  parallel  to  the  coast  and  z  is  positive  upwards.  The 
Navier-Stokes  or  momentum  equations  for  fluid  motion  are  the  governing 
dynamic  equations: 
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where  u,  v  and  w  are  the  components  of  the  water  velocity  in  the  x,  y 
and  z  directions  respectively,   t  is  time,  p   is  the  water  density,  p  is 
pressure,  u   is  molecular  viscosity,  oj  is  the  angular  velocity  of  the 
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Earth,  <^  is  the  latitude,  v^  is  the  Laplacian  operator  v^  =  -— 2- 

a2        32 

+  i~r  +  ^— z  »  3"cl  the  fluid  has  been  assumed  incompressible.     The 

continuity  equation  is: 


|ii+9v^3w  (2-ld) 

ax      ay      az  ^        ' 

Equations   (2-la)  through   (2-lc)  portray  explicitly  only  the  laminar 
processes.     In  this  development,  the  turbulent  effects  are  dominant, 
especially  the  turbulent  or  Reynolds  stresses.     Although  it  is 
impossible  to  describe  the  turbulent  motion  in  detail,  it  is  useful 
to  separate  the  effects  of  the  mean  and  fluctuating  velocity  components. 
By  describing  the  variables  in  Equations   (2-1)  as  the  sum  of  a  time 
mean  value  and  a  fluctuation     and  then  taking  a  time  mean,  the  turbulent 
equations  of  motion  result.     Thus  setting  u  =  Lr+u',v  =  v  +  v', 
w  =  w  +  w'  ,  and  p  =  p"  +  p'    leads  to  the  x-equation  of  motion 


9U^.-3U   +   -au.^.-9U_L32+H.v2    u    -  ( iZI  +  ^V.  ^    aUZZ)     (2-2) 

at    ax    ay    az    p  ax   p     ^ax    ay    az   ^  ^  ^ 
where 


and  Ui,  U2  and  U3  indicate  u',  v'  and  w'  respectively.  The  terms 
in  parenthesis  are  the  turbulent  stresses  and  T  is  a  period  large 
compared  to  the  period  of  a  turbulent  fluctuation  but  small  compared 
to  any  change  in  u   or  7.  It  is  assumed  here  that  terms  of  the  form 
—  v^  IT  are  negligible  when  compared  to  the  turbulent  stresses.  Thus 
the  turbulent  equations  to  be  used  are: 
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!!^"l^^"l^^»lf-M.1n»)v=-l|£.i(^.g^.;^)  (2-3a) 
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where  it  is  now  understood  that  u,  v,  w,  p,  t'    indicate  time  mean 

^i^j 
values,  x'   ,  the  Reynolds  stress,  is  represented  by  -puTuT  where 

pu'-u'.  =pu!.u'..  This  is  explained  in  detail  in  reference  [18].  The 

continuity  equation  in  turbulent  form  is  unchanged 


Frequently  a  Boussinesq  approximation  is  used  for  the  expression 
in  parenthesis  in  Equation  (2-2)  which  maintains  the  original  form  of 
Equation  (2-1).  The  equations  remain  intractable  either  in  the  form 
(2-l)or(2-3).  Since  the  velocity  gradients  in  the  vertical  are  antici- 
pated to  be,  in  most  cases,  an  order  of  magnitude  larger  than  the 
velocity  gradients  in  the  horizontal,  the  lateral  stresses  are 

neglected.  Thus  it  is  assumed  that  x'  =  x'  =  0.  The  normal  stress 

■\^i  xy    yx 

x-x- 
terms  of  the  form   "*  ''  are  assumed  small  compared  to  the  surface 

8xi    1  3p 

pressure  term         -r:;^    and  the  hydrostatic  pressure.  To 

P  oX 

further  simplify  Equations  (2-3)  the  assumption  is  made  that  the 
vertical  accelerations  are  small,  or  as  a  consequence  that  the  pressure 
distribution  in  the  vertical  is  hydrostatic.  Thus, 

p(x,y,z,t)  =  p^(x,y)  +  Y(n(x,y,t)-z)  (2-4) 


n 


where  ri(x,y,t)  is  the  free  surface  displacement  from  mean  sea  level 
(MSL),  p  (x,y)  is  the  pressure  at  the  free  surface,  and  y  =  pg  is  the 
specific  weight  of  water. 

The  Equations  (2-3)  are  now  integrated  over  the  vertical  from 
-h  to  n  where  -h  is  the  z  coordinate  of  the  bottom.  This  procedure 
together  with  Equation  (2-4)  and  the  free  surface  and  bottom  boundary 
conditions  leads  to  the  vertically  integrated  or  tidal  equations  of 
motion.  The  details  carried  out  in  Appendix  A  are  straightforward 
and  with  approximations  to  be  defined  later  lead  to  the  following, 
using  Leibnitz's  rule. 


^%         \    ^%         %    ^%  D    ^Pn  ftn         1 

^'D-^'^W-  2^(--"*)^y  =   -  '  air  -  9°  li  ^  }^\  -b   )    (2-5a) 
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where  the  subscripts,  n,  and  b  indicate  the  subscripted  quantities  are 
to  be  evaluated  at  the  surface  and  bottom  respectively.  D,  q  and  q 
are  defined  as  follows: 


D(x,y)  =  h(x,y)  +  n(x,y,t) 

n(x,y,t) 
qv(x,y,t)  =  I       u(x,y,t)  dz 

'-h(x,y) 

•n(x,y,t) 
qy(x,y,t)  =1       v(x,y,t)  dz 

h(x,y) 


ix(X'y»t)  = 

[x.y.t)  =J 
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It  should  be  noted  also  that  the  convective  terms  of  the  form 

^  T-T-  as  written  are  only  approximate;  the  actual  terms  resulting 

from  the  integration  are  of  the  form 


J  -h 


but  are  written  in  the  differential  form  in  analogy  to  Equations  (2-1). 
The  convective  terms  will  be  neglected  in  the  further  development  of 
the  model  until  considered  again  in  Appendix  B.  In  summary,  the 
equations  of  motion  in  the  final  form  to  be  used  are: 


^  -2»(s1n*)q^  =.£!!".  go  1^  .  1  U^^-.,J  (2-6a) 


The  form  of  the  continuity  equation  is: 


3q    3q 
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b.  Assumptions  and  Boundary  Conditions 

The  equations  of  motion  were  developed  in  the  preceding 
section.  To  specify  a  particular  problem,  various  conditions 
must  be  imposed  on  the  boundaries  of  the  model.  On  the  water 
surface  a  forcing  function  consisting  of  a  normal  pressure  and  a 
shear  stress  is  applied.  At  the  bottom,  a  shear  stress  or  friction 
dependent  on  the  water  velocity  is  imposed.  In  addition  conditions 
on  the  surge  height  and  water  transport  are  necessary  at  the  areal 
boundaries  of  the  model. 

The  shear  stress  at  the  free  surface  is  calculated  by  the 
method  proposed  by  Van  Dorn  [19].  The  shear  stress  is  based  on  the 
wind  velocity  and  is  given  as  follows: 


^.    =  pMu|u, 


x,^=pk|U|U^ 


(2-7) 


whe 


re  U  =    /u^  "*■  U,^  and  U     and  U     are  the  x  and  y  components  of  the 
V    X  y  X  y 

wind  velocity,  respectively,  at  a  reference  elevation  of  approximately 
10  meters.     The  constant  k  has  been  determined  by  Van  Dorn  as  the 
following  function  of  velocity: 


;^    ;  U  <     U 


cr 


k  =/ 


(2-8) 


ici   ^  4(1    -  U^^/U)^    ;  U  ^  U^^ 
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where 

k^  =  1.1  X  10"S 

k^   =  2.5  X  10"^ 

U^^  =  14  knots  =  23.6  feet/second 
cr 

There  are  uncertainties  in  applying  Van  Dorn's  results  to  hurricane 
winds.  Van  Dorn's  data  are  generally  for  wind  speeds  low  in  comparison 
to  the  wind  speeds  of  a  large  hurricane.  Although  there  is  no  proof 
that  the  results  are  directly  applicable,  they  are  used  here  because 
the  associated  measurements  showed  little  scatter  and  the  results 
are  reasonably  consistent  with  other  data  sources.  VJu  [20]  has  pub- 
lished more  recent  information  concerning  the  wind  stress  coefficient, 
k.  Figure  2-2  summarizes  available  data.  Van  Dorn's  k,  calculated 
from  Equation  (2-8),  is  also  plotted  in  Figure  2-2  to  provide  a  com- 
parison. Since  the  hurricane  model  to  be  discussed  requires  a  shear 
stress  field  as  input,  a  better  shear  stress  model,  if  it  becomes 
available,  can  be  easily  incorporated  into  the  surge  model  without 
changes  in  concept  or  approach.  The  wind  velocity  input  will  be 
discussed  in  Section  2-c. 

The  treatment  of  the  bottom  friction,  f.  ,  t,   is  a  much 

X   y 
more  difficult  problem.  In  making  storm  surge  calculations  the 

simplest  assumption  is  to  consider  x,   =  t.  =  o.  This  is  employed 

X    y 
in  many  developments  of  storm  surge  calculations,  especially  the 

earlier  ones.  Under  appropriate  conditions  this  is  a  justifiable 
assumption.  For  the  case  of  a  rapidly  moving  storm  traveling  per- 
pendicular or  nearly  perpendicular  to  the  shore  it  is  a  good  assump- 
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tion  since  there  is  little  time  for  "wind-up"  and  the  friction 
effects  are  small  due  to  low  water  velocities  [13]. 

A  better  assumption  in  some  cases  is  that  of  a  linear  or 
"linearized"  bottom  stress  of  the  form  fiq  ,  f.q  where  f,  is  a 
linearized  friction  factor.  This  is  the  type  of  friction  used,  for 
example,  by  Dean  and  Pearce  [21].  A  "linearized"  friction  term  is 
generally  necessary  if  friction  is  required  and  any  type  of  analytical 
solution  is  desired.  The  bottom  shear  stress  employed  in  this  model 
is  the  Darcy-Heisbach  expression  and  is  the  most  accurate  of  the 
three, 

P^hlu^  _  pf|u|Uy 

^b^  "    8    '  ^b  '    8 

X  V 


or  in  terms  of  volume  transport 


pflqhx  pf|q|qy 

^b  "   8D^    '  ^b„  "   SD-^ 
A  y 


where  u  =  /u  ^  +  u  -  ,  q  =  /q  ^  +  q  ^  ,  and  f  is  the  quadratic 
bottom  friction  factor.  In  incorporating  this  steady  state  rela- 
tionship into  the  model  it  is  assumed  that  the  law  is  valid  for  the 
propagation  of  long  waves  of  the  period  of  interest.  This  is 
equivalent  to  assuming  that  the  velocity  variation  over  depth  is 
nearly  the  same  for  slowly  varying  flows  as  for  steady  state. 
Although  this  technique  has  limitations,  especially  in  deep  water, 
it  is  the  best  method  available  that  can  be  used  in  conjunction  with 
the  nonlinear  equations  to  allow  for  flooding  and  to  permit  cal- 
culations in  very  shallow  water.  Another  important  consideration  is 
that  in  deep  water  the  friction  is  relatively  less  important  than  in 


17 


shallow  water  since  the  water  velocities  are  in  general  much  smaller. 
As  an  illustrative  example,  the  response  of  a  water  column  to  an 
impulsive  wind  is  depicted  in  Figure  2-3. 

Returning  to  the  problem  specification,  the  boundary  conditions 
are  illustrated  in  Figure  3-5.  At  the  seaward  boundary,  the  wind 
setup  is  assumed  zero  as  a  result  of  considering  an  infinite  depth 
at  that  point.  The  barometric  tide  still  exists  at  the  seaward 
boundary  and  is  given  by  n(x,y)  =  1.15(p^-p  (x,y))  where  p^  designates 
the  barometric  pressure  at  a  distance  large  compared  to  the  storm 
dimension,  where  n  is  in  feet,  and  the  pressure  p  is  in  inches  of 
mercury. 

At  any  land-sea  boundary  the  flux  normal  to  that  boundary  is 
zero,  q^  =  0,  where  the  subscript  n  indicates  the  component  normal 
to  that  boundary. 

Most  difficult  from  a  physical  standpoint  are  the  lateral 
boundary  conditions.  They  are,  however,  necessary  to  carry  out  the 
calculations.  Since  no  real  physical  limitations  exist,  conditions 
must  be  chosen  that  least  disturb  the  model  and  therefore  allow  a 
realistic  solution  to  be  obtained.  A  further  condition  is  then 
imposed  that  the  boundaries  are  at  a  distance  large  compared  to  a 
characteristic  dimension  of  the  storm.  This  means,  in  effect,  that 
the  boundary  conditions  on  the  lateral  boundaries  become  a  second 
order  effect  if  they  are  far  away.  The  boundary  conditions  used  at 
the  lateral  boundaries  of  this  model  are:  the  flow  in  the  x  direction 
is  equal  to  zero,  q  =  0,  and  the  surge  height  n  is  set  equal  to  the 
barometric  tide  at  the  lateral  boundaries. 
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The  initial  conditions  are  simply  a  quiescent  system.  The  surge 
heights  can  be  initially  zero  or  in  barometric  equilibrium;  this  makes 
little  difference  since  the  storm  is  initially  far  from  the  region 
of  interest.  Initial  barometric  equilibrium  is  used  in  this  study. 
More  importantly  the  storm  is  assumed  to  start  at  a  distance  offshore 
large  enough  that  the  effects  due  to  model  start-up  are  negligible. 
The  effects  of  various  starting  points  can  be  seen  in  Figures  4-5  and 
4-10,  where  the  surge  height  is  relatively  insensitive  but  the  currents 
exhibit  wind-up. 
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c.  The  Hurricane  Model 

The  actual  storm  surge  model  uses  t  ,  t   and  p  as  input. 

Since  the  measured  wind  shear  stress  is  in  general  not  available, 

Van  Dorn's  technique  is  used  to  convert  U  and  U  to  t   and  x   as 

X     y    n^     ny 

discussed  in  Section  2-b.  Actual  wind  fields  may  be  applied  or  an 
idealized  model  used  to  represent  the  storm.  In  this  section  an 
analytical  model  is  used  to  represent  wind  velocity  and  pressure  fields. 
Since  the  actual  surge  model  is  the  primary  objective  of  this  study, 
the  hurricane  model  will  be  presented  and  used  without  further  comment, 
except  for  a  short  note  on  the  numerical  technique  in  Chapter  3. 

The  hurricane  representation  used  is  not  unique  and  a  different 
one  could  be  substituted  into  the  numerical  model.  The  hurricane 
model  used  is  that  presented  by  Wilson  [23].  The  important  parameters 
of  an  idealized  hurricane  are:  the  central  pressure  index  CPI  or 
equivalently  the  central  pressure  anomaly  Ap  =  p_^  -  p  (where  p^  is 
the  normal  pressure  or  the  pressure  at  a  large  distance  from  the 
hurricane  center  and  p  is  the  central  pressure  of  the  hurricane),  the 
distance  from  the  hurricane  center  to  the  point  of  maximum  winds  R, 
and  the  forward  velocity  V  of  the  hurricane.  These  parameters  Ap, 
R,  and  V  are  reported  by  the  National  Weather  Service  when  character- 
izing a  storm,  thereby  providing  a  method  to  approximate  a  hurricane. 

The  isobars  are  assumed  circular  and  the  pressure  distribution  at 

-R/r 
the  surface  is  taken  as  p  =  p^  +  Ape  '  .  The  streamlines  are 
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assumed  circular  for  a  stationary  storm.     The  equations  presented 
include  the  effect  of  a  constant  translational    velocity,  V, 
superimposed  on  the  storm  as  shown   in   Figure  2-4.     Using  Wilson's 
notation: 


where 


c         q 


and 


Up  =  /Ap.    R     -R/r 


V 


c     ./  ^^    -  e 
a 

Ap     R  „-R/r 

U     =  ^°,    , 
g      2w(sin())) 


V   =  V(sine) 


The  subscript  G  refers  to  gradient  wind,  c  to  cyclostrophic 
wind,  and  g  to  geostrophic  wind,  p^  is  the  density  of  air,  V  the 
forward  velocity  of  the  storm,  e  is  defined  in  Figure  2-4  and  2co(sin<)) 
is  the  Coriolis  parameter  as  used  in  Equations  (2-1).  The  gradient 
wind  Ug  is  the  desired  parameter.  These  winds  are  calculated  at  a 
point  high  enough  above  the  surface  to  allow  the  assumption  of  zero 
friction.  To  reduce  these  to  the  10  meter  elevation  a  reduction  factor 
is  needed.  The  factor  depends  on  latitude,  divergence  of  wind 
direction  from  the  isobars  and  other  factors.  In  this  study  a 
constant  factor  of  0.83  was  used.  See  references  [24,25]. 
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Up  =  Gradient  Wind  Velocity 


V  =  Forward  Velocity 
of  Storm 


FIGURE  2-4  DEFINITION  SKETCH  FOR  STORM  TRAVELING  AT  CONSTANT  VELOCITY 


The  winds  of  Hurricane  Camille  are  evaluated  using  this  method. 
The  hurricane  position  as  a  function  of  time  is  given  in  Figure  2-5. 
The  central  pressure  p  ,  and  the  peripheral  pressure  p^,  are  given  as 
functions  of  time  in  Figure  2-6a.  Figure  2-6b  depicts  radius  to 
maximum  winds  R,  and  fonvard  speed  of  the  storm  V,  as  functions  of 
time.  Figure  2-6c  gives  maximum  sustained  wind  speeds  in  knots  as 
a  function  of  time.  These  data  are  from  the  National  Weather  Service 
[26,27].  Figures  2-7a  through  2-7c  present  comparisons,  at  different 
times,  of  the  surface  wind  fields  for  Hurricane  Camille  as  developed 
by  the  National  Weather  Service  [27]  with  the  wind  fields  calculated 
by  the  above  model . 
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7th 


FIGURE  2-5    STORM     TRACK    FOR    HURRICANE 
CAMILLE,     AUGUST    16-18,     1969 
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FIGURE  2-6a 


PRESSURE      CHARACTERISTICS      OF 
HURRICANE      CAMILLE      AS      A 
FUNCTION       OF       TIME 
(DATA      FROM      REFERENCE     [26]) 
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FIGURE  2-6b      RADIUS      TO      MAXIMUM      WINDS      AND 
FORWARD    SPEED      OF     HURRICANE 
CAMILLE      AS      FUNCTIONS     OF      TIME 
(DATA      FROM      REFERENCE     [26]) 
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FIGURE  2-6c 


MAXIMUM      WIND     SPEED     IN     HURRICANE 
CAMILLE      AS     A     FUNCTION      OF     TIME 
(DATA      FROM      REFERENCE     [26]) 
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d.  Analytical  Treatment-Steady  State,  Without  Bottom  Friction 

i .  Problem  Definition 

Having  defined  the  problem  and  established  equations  of  motion 
and  boundary  conditions,  it  remains  to  find  a  solution  for  a  particular 
storm  and  affected  coastal  region.  Before  presenting  the  numerical 
approach  necessary  to  attain  a  solution  of  complicated  prototype  con- 
ditions and  geometry,  several  simplified  analytical  approaches  will  be 
considered.  These  will  serve  three  basic  purposes.  First,  they 
allow  highly  simplified  calculations  for  basic  engineering  purposes. 
Second,  they  allow  intuitive  investigations  into  the  physical  processes 
involved.  Third,  they  produce  solutions  with  which  to  check  the 
numerical  approaches. 

Figure  2-8  is  a  definition  sketch  for  this  section.  The  present 
treatment  includes  only  motions  in  the  x-direction;  Equations  (2-5a) 
and  (2-5c)  therefore  reduce  to  Equations  (2-9a)  and  (2-9b) 
respectively. 


W^  D--X ^3^-9D^V  ^  V  ■  'b  )  ^2-9a) 

A        A 

3q    . 
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FIGURE  2-8     NOTATION   FOR  ONE-DIMENSIONAL  PROBLEM 


The  solution  is  novi  independent  of  y.     Two  cases  will   be  considered, 

the  uniform  depth  case  h(x)  =  h  =  constant  and  the  case  of  h(x)  =  h     - 

mx  or  uniform  slope.      In  both  cases  it  will   be  assumed  that  the  wind 

stress  on  the  surface  is  uniform.      For  the  steady  state  case,   the 

case  of  no  net  transport  in  the  x-direction,  q^^  =  0  and  therefore  the 

shear  stress  t,      =  0,   regardless  of  f,   the  bottom  friction  coefficient. 

x 
(Note:     it  is  assumed  here  that  there  is  zero  velocity  over  the  entire 

q.  s^x 

depth.)   For  this  problem  the  p—  {- — )  term  is  zero  and  the  atmospheric 

3p         °   ^^ 
pressure  gradient  — —  is  assumed  small  in  relation  to  the  other  terms. 

oX 


ii .  Uniform  wind  stress  Acting  Over  a  Shelf  of  Uniform  Depth 

In  the  first  example  h(x)  is  uniform  (=h)  as  shown  in  Figure 
2-9. 
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FIGURE  2-9  ONE-DIMENSIONAL  SHELF  WITH  UNIFORI'l  DEPTH 


The  further  assumption  is  made  that  the  shear  stress  is  uniform  and 
constant  at  the  water  surface  z  =  n.  The  problem  can  be  stated  as 
follows: 


^q^ 


qv  = 


X   at 


Since  D  =  h  + 


pg(h.n)  IJ  ==  x^ 


(2-10) 


Rewriting  Equation  (2-10)  leads  to 


an  _  "x 
9x  ~  pg(h+n) 


(2-11) 


At  X  -  0,  the  boundary  condition  for  n  is  determined  by  noting  that 
seaward  of  that  point  h  ->  »  and  therefore 
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1^-^  0    as     h  -V  CO 

Thus  n  =  constant  as  h  ^  «>,  and  n  =  0  is  chosen  as  the  boundary  con- 
dition at  the  shelf  break,  x  =  0.     Solving  for  n  yields 

T 

n(h  +  n/2)  =  Ax 


or 


I        \ 


n  =  -h  +    Vh2  +  _^2x 

pg 


Evaluating  x     ,  using  Van  Dorn's  constants,   completes  the  solution. 

A 

An  even  simpler  procedure  is  frequently  the  basis  for  storm  surge 
calculations.  Assuming  that  the  storm  surge,  n,  is  small  compared  to 
the  undisturbed  depth  (i.e.  h  +  n  "^  h)  leads  to  the  solution 

T 

n  =  -rf  X 


P 


W 


Designating  x  as  the  fetch  or  shelf  width  and  h  as  depth  leads  to 
the  storm  surge  at  the  coastline  x  =  x 

T 

I       ^x   Fetch 
n 


x=x    pg    Depth 

Figure  2-10  provides  a  simplified  procedure  for  making  these 
calculations  and  is  based  on  the  latter  solution.  As  an  example, 
assume  a  depth  of  50  feet  for  a  50  mile  shelf  width  (or  fetch)  and 
a  wind  speed  of  50  knots,  then  the  storm  surge  calculated  from 
Figure  2-10  is  3.23  feet. 
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1 i i .  Uniform  Wind  Stress  Acting  Over  a  Shelf  of  Uniform  Slope 

As  a  slightly  more  realistic  situation,  the  same  assumptions 
are  made  as  in  the  previous  case  except  now  a  sloping  bottom 
(or  h(x)  =  h  -  mx)  is  considered  as  shown  in  Figure  2-11. 


♦  X 


FIGURE  2-11  ONE-DIMENSIONAL  SHELF  WITH  UNIFORMLY  SLOPING  BOTTOM 


From  Equation  (2-11 ) 


(h  +  n)  Iv  =  -^  =  constant  =  e 
dx   pg 


and  noting  that  h(x)  =  h  -  mx 


D  |2  =  (h.n)  If  .  (ht„)  1^  .  (h.,)  |a  -  „  (h.,) 


a.u   J    DdD 

at  X  =  0  to  obtain 


Integrate  and  apply  the  boundary  condition  n  =  0 


1  -  C, 


4  =  in  (rr— -,) 


(2-12) 
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where 


r  =  OH 


mh 


mh 
mhi 


To  simplify  the  calculations  obtain  n  at  the  coastline  x  =  x  . 


Equation  (2-12)  then  reduces  to 


E  =  1"  (1  !  '°  ■  ;)  C2-13) 


Solutions  to  Equation  (2-13)  at  x  =  x  are  provided  in  Figure  2-12. 
As  an  example  consider  a  shelf  with  x  =  50  nautical  miles,  h  =  80 
feet,  and  h,  =  20  feet,  this  shelf  has  the  same  width  and  the  same 
average  depth  as  in  case  ii.  For  the  same  50  knot  wind,  n  at  the  coast 
is  3.53  feet.  Heuristically ,  one  would  expect  a  larger  n  for  the 
case  of  the  sloping  shelf. 

Consider  the  location  on  the  sloping  shelf  which  has  the  same 
depth  as  the  uniform  depth  shelf.  For  a  point  Ax  towards  shore 
there  is  a  change  in  depth  Ah  and  conversely  for  Ax  away  from  shore. 
Thus 

^^+  ^  constant  ^  constant  ^  .  Ah  +  (Ah^^  _  ^  ^^^^3^  . 

An"  =  constant  ^  constant  m  +  Ahi  +  (Ahj^  ^  ^  ^^^^^3^ 
^^^  K  FT    IT 
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and  the  total  over  2ax  is 

6n  =  An^  +  An"  =  S^nilM.  [2  +  2(^)2  +  0  (Ah)'^] 

It  can  be  seen,  then,  that  a  shelf  of  uniform  depth  will,  in  all 
cases,  produce  a  smaller  surge  than  a  shelf  of  constant  slope  with 
the  same  average  depth. 
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e.  Analytical  Treatment  -  Moving  Wind  Stress  and  Pressure  Field 

i.  Problem  Formulation 

The  purpose  of  this  section  is  to  contribute  to  the  development  of 
a  simple  procedure  for  the  analysis  of  the  dynamic  effects  involved 
in  storm  surge  calculations.  This  analysis  uses  the  one-dimensional 
tidal  equations  and  the  development  will  be  for  an  arbitrary  forcing 
function  such  as  a  hurricane.  Also  the  bottom  profile  must  be  in  an 
analytical  form  simple  enough  to  allow  a  solution  to  the  equations  of 
motion.  Conceptually,  the  approach  is  as  follows:  find  the  storm 
surge  ri(x,t)  resulting  from  an  impulsive  forcing  function  representing 
the  wind  shear  stress,  pressure  gradient  or  both.  If  the  impulse 
occurred  at  time  t'  and  position  x',  the  resulting  surge  at  (x,t) 
could  be  written  n(x,t)  =  G(x,t;  x',t');  this  is  sometimes  referred 
to  as  the  Green's  function.  Having  found  the  general  form  of  the 
solution  G(x,t;  x',t')  for  a  simple  offshore  bathymetry  one  can  then 
construct  the  solution  for  an  arbitrary  forcing  function.  This  may 
be  obtained  by  formal  integration  or  the  results  may  be  obtained 
numerically. 

The  development  in  this  saction  v/ill  be  for  the  special  case  of 
a  shelf  of  finite  width  and  uniform  depth  as  shown  in  Figure  2-13 
and  a  triangular  shear  stress  distribution  as  shown  in  Figure  2-15  . 

For  the  equations  of  motion  and  continuity,  consider  Equations 
(2-5)  where  q  is  abbreviated  as  q,  t   as  t  ,  and  x^^  as  t^^. 


39 


l^^&lf^9°lf=-^^^^(^-b)    (2-14a) 


||-+|a=0  (2-14b) 

The  convective  acceleration  terms  will  be  neglected  and  the  assumption 
n  <  <  h  is  made,  and  since  D(x,t)  =  h  +  n(x,t)  ^   h,  Equations  (2-14a) 
and  (2-1 4b)  reduce  to 

9q(x.t)   9n(x,t)  ^  n 

3x     n 

Ignoring  the  bottom  friction  term  x,    and  letting   (-  x   (x,t)   - 

h   ^Pn(^'t)  '     ' 

-r^^^ )  =  F'(x.t)  we  have: 

|a+gh  |i=  F'{x,t)  (2-15a) 

fa +1^=0  (2-15b) 

The  variables  are  described  in  Figure  2-13,  where  q  is  the  volume 

transport  per  unit  width.     Differentiating  Equation  (2-15a)  with 

respect  to  x  and  Equation   (2-1 5b)  with  respect  to  t  and  subtracting, 
produces  the  inhomogeneous  wave  equation 

^^^  „k       9^1^      -  8F'(X,t)  r-f  .X  .  fr,       TC\ 
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FIGURE  2-13  IDEALIZED  SHELF  BATHYMETRY 


i  i .  Solution 

First  a  solution  to  the  homogeneous  problem  F(x,t)  =  0  is 
obtained.  Using  separation  of  variables,  let  n  =  X(x)T(t) 

For  the  spacial  part  of  the  problem,  we  obtain 


-gh 


3F-  -  V' 


(2-17) 


where  a  '^   is  a  constant  as  yet  undetermined.  The  solution  to 
Equation  (2-17)  is 


a  a 

Ci  sin  X  +  C2  cos  — ^  x 


(2-18) 


/gH"    ^   /gh 

Referring  to  Figure  2-13  note  that  the  depth  is  infinite  at  x  =  0"^ 
and  the  surge  must  be  zero  at  that  point.  At  x  =  z,   there  can  be  no 
flow  through  the  wall  and  this  leads  to  the  condition  ^  =  0  at 
X  =  £. 
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The  condition  n  =  0  at  x  =  0  determines 

C2  E  0  and  1^  I    =  — ^  Ccos  -^)   T(t)  =  0 

a  £ 
leading  to  - —  =  (2n  +  1)  5-  .  The  general  case  for  the  solution 

becomes  a  sum  of  the  solutions  in  Equation  (2-18).  Thus  the  solution 

is  taken  to  be  in  the  form 


n  =    E    a^(t)  sin  g^  (2-19) 

n=l,3,5'- 

Referring  to  the  homogeneous  form  of  Equation  (2-16),  obtain  with 
Equation  (2-19)  the  following: 

"   a^a  (t)  9 

2   ^sF ^^"Zr-^^  (21)  ^n^^^  '^"  2^^  -  ° 

n=l ,3,5  ■•  • 

From  orthogonality,  obtain 


a„"(t)  +  0  2  a  (t)  =  0 
n       n   n 


where  a  =  S^^  and  c  =  /  gh 


The  two  homogeneous  solutions  are  of  the  form  a_(t)  =  sin  at, 

00 

cos  ot.     Writing  F(x,t)  =    l         B„(t)  sin  ^  and  B  (t)  = 

2  r^        niix        n=l,3,5  "       ^*-      " 

-  F(x,t)  sin  21^  "^^  ^'^'^  '^^^   "»  ="^^  using  the  method  of  variation  of 

parameters,  arrive  at  the  following  for  the  inhomogeneous  equation; 

with  the  lower  limit  determined  by  the  fact  that  a  (0)  =  0. 

,    X       ^^"  ""n^        (   ^  coso^t  p 

^n^^)  =         a„  ]         B^(t')  cosa^fdf    --^        B,^(t')sina^t'dt' 
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a,(t)  = 


Bjf  ) 


sina^(t-t')dt' 


(  t 


E     a^(t)  sin  ^ 


mix 


n=l  ,3,5 


n=l  ,3,5  °n  j 


B^(t')sina^(t-t' )dt' sin  g^ 


and  finally 


n(x,t)  =   E 


n=l,3,5  % 


%   ft 


0  J 


F(x',t' )sina^Ct-t')sin  ^^  dt'dx'sin^  (2-20) 


Letting  F(x',t')  =  k6(x'-x  )  6(t'-t  )  one  then  has  an  impulse  function 
located  at  x  =  x  and  t  =  t  where  n  represents  the  Green's  function 
or  response  to  this  impulse  function,  and  <   is  a  constant  representing 
the  strength  of  the  instantaneous  surface  shear  stress  function 


n(x,t)  =  G(x,t;  x^.t^)  = 


n=l,3,5  ^''n 


mrx 


k6(x'-Xq)  6(t'-tQ)  sina^(t-t')sin  ^— -dt'dx'sing- 


nTrx 


OJ 


niTX. 


G(x,t;x^,t^)  =   ^  J^  sina^  (t-t^)  sin  -^  sin  g^ 
n= 1,3,5   n 


(2-21) 


It  is  now  possible  to  construct  a  solution  by  the  superposition  of  an 
appropriate  series  of  impulses  in  space  and  time. 

For  some  hurricane  motions,  a  more  convenient  form  of  the  Green's 
function  is  obtained  by  using  an  impulse  function  which  travels  with 
constant  velocity  and  remains  unchanged  during  its  traverse  over  the 
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shelf.  Assuming  the  forcing  function  to  always  start  at  x  =  0,  to 
travel  at  constant  velocity  C,  and  that  the  time  that  it  crosses 
X  =  0  is  arbitrary  at  t  ,  the  impulse  function  is: 


F(x,t)  =  k6(x  -  CCt-t^))  0  £  X  £  Jl 


Without  loss  of  generality  let  t  =  t  ^  t  which  yields 


F(x,t)  =  k6(x  -  Ct) 


0  <  X  <  £ 


Thus  for  the  case  of  a  suddenly  increased  shear  stress  as  shown  in 

T 

Figure  2-14,  <  =   ,  where  t  is  the  amount  of  the  increase  in 

shear  stress. 


T  S(x-Ct)  =  F(x,  t) 
P 


-^C 


-►  X 


FIGURE  2-14  IMPULSIVE  FORCING  FUNCTION  MOVING  AT  CONSTANT  SPEED 
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For  the  case  of  an  arbitrary  shear  stress  distribution  as  shown  in 
Figure  2-15a,  locally  Ax  =  tt  Ax 

oX 


♦-X 


FIGURE  2-1 5a     ARBITRARY   FORCING   FUNCTION 


For  a  triangular  shear  stress  distribution  as  shown  in   Figure  2-15b 


-Kq 


>   X 


♦•  X 


FIGURE  2-1 5b  TRIANGULAR  FORCING  FUNCTION 
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K  becomes  —  tt^  5x,  and  the  following  is  obtained. 

p    3X  ^ 


F(x',t')    =  i|j  6(X'-Ct')    dx    =  -  At    5(X'-Ct') 


The  resulting  Green's  Function  can  then  be  written  in  integral   form  as 
follows. 


Cl  Ct 


G(x.t;0,0)=       E       -^ 
n=l  ,3,5     n J 


At    ,/..,    ^j.,>_.-^       rj.    ^1  \_.--.nTTX' j^i  j..,^  .    nirX 


6(x'-Ct')sina„(t-t')sini^^^t'dx'sin^ 


OJQ 


Note  that  this  is  really  the  solution  for  an  incremental   element  having 
a  shear  stress  jump  of  At.   Since  6(x'-Ct')  f  0  only  when  x'   =  Ct' 


G(x,t;0) 


2At 
n=lt3.5  P% 


/j.    j.i\     •    nirCt'      ,.|       .       nirX 

sina^(t-t  )sTn-2^^ — dt'    sin  2;^ 


■     ...  nirC  I  I         nirC 

Letting  a^   =  a  =  ^  ,  o^     =  a     =^ 


and  integrating  yields 


2a 


.•-  nTTX 


G(x,t;0)  =       z  ■    /    ,  2 — ZT    ^^'    ^■'"'^'t  "  °  sina't]  sin  ^ 

n=l,3,5  P^°^°     "°   '' 

Note  that  this  is  the  solution  for  a  solitary  continuous   impulse 
function  traveling  at  speed  C;  at  time  t  >  z/C  the  forcing  function 
crosses  the  shelf,  and  the  solution  becomes  periodic  in  time. 
Letting  c*  =  C/c,  c  =    /gh  and  a'   =  a  c*  one  arrives  at 


+2at 


-•-  n-iTX 


G(x,t;0)  =      E        pzoHc*^-])   [c*sinat-sin(ac*t)]  sin  ^ 


n=l  ,3,5 


with  the  constraints  t  <  z/C  and  c*  M  •     When  t  >   ii/C 
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G(x,t;0)  =    z 


2ai 
pZa 


rz/C 


n=1.3,5  P^°  J 


sina(t-t')  sin  ^df  sin  g^ 


where  a'   i/C  =  ^ 


and 


COS  2~  ~  0 

n-T 
sin  ^=   (-1)  2 


n  =  1,3,5 


G(x,t;  0)  thus  becomes 


n  T 

Gfx  t-f)^  =         y         2ATSin(nTix/2£)r,    ^   .  ,  -5— 

n=l'3,5  ^I^^#^^1T~^'^^^  ^"""^t  -  cosa(t-ii/C)    (-1)2] 

with  the  constraints  t  >  £/C  and  c*2  ^  1 .     if  at  is  written  as  ^  dx' 
and  it  is  noted  that  x  =  Ct  then  At  =  C  |j  dt.     Letting  t  =  t  -  t 

dX  Q 

yields  G(x,t;  t^)  where  t  represents  the  time  at  which  6(x-C(t-t  )) 

0 
crosses  the  point  x  =  0.  The  desired  result  is  then  the  sum  or  integral 

of  the  above  solutions 


n(x,t)  =  1    G(x,t;  t^)  dt^ 
0 


(2-22) 


The  shear  stress  profile  of  a  hurricane  can  be  approximated  by   a 
triangular  function  [2].  The  derivative  of  this  is  two  rectangular 
pulses  as  indicated  in  Figures  2-15b.  The  response  can  be  found 
using  Equation  (2-22)  and  remembering  that  G(x,t;  t^)  is  represented 
by  two  functions  depending  on  whether  t  is  less  than  or  greater  than 
i/C.     The  solution  is  straightforward  but  unfortunately  tedious. 
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The  solution  is  simplified  by  putting  n  into  a  dimensionless  form. 
Letting 

T  =  4JI/C 
a*  =  aT  =  2niT 
X*  =  x/£ 
t*  =  t/T 

allows  n'   to  be  represented  as  Equation  (2-23). 


•rv*  i-*\  =  m 


n'Cx*,t*) 


(X,t)    (2Tr)3 


X27TVI7 


(2-23) 


The  appropriate  variables  and  constants  for  the  triangular  stress 
distribution  are  illustrated  in  Figure  2-15,  where  ti   has  been 
arbitrarily  set  equal  to  zero  for  convenience. 
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FIGURE  2-16  DEFINITION  SKETCH  FOR  SHEAR  STRESS  DISTRIBUTION 


and  are  defined  as  follows 


k'  =  -  — 

p  3X 
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ti  =  0  -V  ti*  =  0 


t2    = 

X 
2C 

>   t2* 

X* 
8c* 

t3    = 

X 

t3*   = 

X* 

4c* 

X*  = 

x_ 

K*  =  -<;-,   0  <  t^  <  1^ 


I  _   I      X      1      X 

0    2C    0   C 


..  _  o  .   X 


C    0 

The  solution  must  be  presented  as  six  different  forms,  each  for  a 
different  time  span.  They  are  as  follows: 
(1)  In  the  first  case  the  storm  has  completely  passed  by  the 
position,  X  =  Ji,  thus 

.  x+i       .J.   x*+l 

and  the  solution  for  n  is 

n'(x*,t*)  = 

n-1  tQ=t3 

J:    'jfte'(f-^  [sina((t-tJ-ii/C)  (-1)  2  +c*cosa  ft-t  )] 
n=l,3,5  "  ^c  -'^       n    0  "    °    . 

^o"^l 


By  defining  the  contents  of  the  brackets  as  H-jCt^),  the  response  becomes 
n'(x*,t*)  =  _T.         ^.^[5^[H^(t3)  +  H^Ct^)  -  2H^(t2)] 
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where  t-j ,  t2,  t^  represent  values  of  t  and  are  defined  in  Figure  2-16. 
(2)  In  case  2  the  storm  is  completely  on  the  shelf  0  <  x  <  £,  thus 

X/C  <  t  <  £/C  ^  |^<  t*  <  ^ 

and  n  can  be  written  as  follows 

n'(x*.t*)  = 

*o"^3 

„_,%  ,  n'(^"'-lf'  [^c*cos(o„(t-t„))-  1^  cos(c*a„(t-t„))] 
n  I , J ,D  .  _ 

Defining  H^Ct  )  as  the  contents  of  the  brackets,  as  above,  leads  to 

n= 1 ,3,5        ' 

The  time  periods  from  0  <  t  <  A/C  and  i/C  <  t  <  -x-  have  not  yet  been 

covered  since  they  introduce  the  additional    complication  of  having 

only  part  of  the  forcing  function  in  the  region  of  interest.     The 

solutions  for  these  remaining  time  periods  are: 

(3) 

£  +  A/2       .        a  +  \         1  +  x*/2       .*       1   +  A* 

n'(x*,t*)  = 

2         nJ(c"Mf'^   [H2(t3)-H2(t-£/C)+H^  (t-il/C)-2H,  Ct2)+H^  (t^ )] 

(4) 

./f^   ^  ^.    ^   Jl   +   A/2  1  .^        1    +A*/2 
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n'(x*.t*)  = 


E     ^^3(^1^.(f  ^  [H2(t3)-2H2(t2)  +  H2(t-il/C)-H^  (t-^/C)+H^  (t^ )] 
n=l ,3,5, . . .  ^ 


(5) 

A      1     A     A       J.  ^     A 

n'(x*,t*)  =  ^^^  3^5    ^T(^5?l-y--  [H2(t)+H2(t^)-2H2(t2)] 


(6) 

X* 


0  <  t  <  X/2C  ^  0  <  t 


-  "  -  8c^ 


II      I    }«J}S/)**> 


n'(x*,t*)  has  now  been  described  for  all  t*.   It  can  be  seen  that  n* 

can  be  described  as  a  function  of  three  dimensionless  quantities 

c*.  A*,  and  t*. 

Figure  2-17  shows  profiles  of  n'  on  the  shelf  for  various  values 

of  t*  and  particular  values  of  X*  and  c*. 

Figure  2-18,  perhaps  of  more  interest,  shows  n'   ,  the 
^  max 

maximum  storm  surge  that  occurs  as  a  result  of  the  storm  passage  as  a 
function  of  both  c*  and  x*. 

To  determine  the  accuracy  of  both  theory  and  calculations,  an 
independent  approach  was  utilized  to  obtain  equivalent  dimensionless 
surge  heights.  The  one-dimensional  tidal  equation  of  motion  and 
equation  of  continuity  in  finite  difference  form,  as  developed  in 
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X' 
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FIGURE     2-18 


DIMENSIONLESS  SURGE     HEIGHT 

MAXIMA     AS     A  FUNCTION     OF 

DIMENSIONLESS  STORM     WIDTH     AND 
VELOCITY 
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Section  3-a  ,  were  used  to  obtain  a  numerical  solution  to  this  problem. 
The  equations  were  linearized  (D  =  h)  to  agree  with  the  assumptions 
in  this  section,  and  the  bottom  friction  neglected.  The  "no  flow" 
boundary  condition  q  =  0  was  prescribed  at  the  coast  with  a  boundary 
condition  n  =  0  assumed  at  the  shelf  break.  Note  that  the  q  =  0 

A 

condition  is  comparable  to  ^  =  0.      Figure  2-19a  depicts  a  shelf  of 
constant  depth  divided  into  20  one-dimensional   grid  elements,  as   used 
for  the  calculations.     The  "storm"  or  triangular  shear  stress  dis- 
tribution as  illustrated  in  Figure  2-19b  was  moved  over  the  shelf 
waters  at  a  constant  speed,  C.     A  slight  approximation  was  used  here. 
Once  the  edge  of  a  shear  stress  step  encounters  a  new  grid  element, 
the  stress  was  assumed  constant  over  that  element  for  a  time  At  =  ax/C, 
or  the  time  for  the  leading  edge  to  traverse  the  grid  element.     The 
numerical   results  were  put  into  dimensionless  form  as  in  Equation   (2-23) 
and  compared  to  the  analytical    results  yielding  Figure  2-20. 
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Notes: 

(I)     Tj    Linearly    Extrapoloted     to     Coost 


(2)    Boundary     Conditions 
(o)    q  =  0  at    Coast 

(b)    t;  =  0  in    Grid    (I) 


Shear     Stress    Distribution 
Moves     at      Constant 
Velocity     C. 


20    19      18     17     16      19     14     13      12     II      10      9      8      7      6      S       4     3      2    i  =  l 


•20    Grid     Elements 
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FIGURE  2-l9a    DEFINITION     SKETCH     FOR     NUMERICAL 
CALCULATIONS 
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FIGURE   2-l9b     FINITE     DIFFERENCE    MODEL   OF 

TRIANGULAR    SHEAR    STRESS    DISTRIBUTION 
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.3 
1 


.4      .5    .6    .7  .8  .9  1.0 
(Analytic) 


FIGURE  2-20     COMPARISON     OF     NUMERICAL    AND 

ANALYTICAL  SOLUTIONS 


CHAPTER  3 
THE  NUMERICAL  MODEL 
a.     One-dimensional   Numerical   Model 

In  most  prototype  situations  it  is  unrealistic  to  make  the  simpli- 
fying assumptions  necessary  for  a  successful  analytical  treatment  of  the 
storm  surge  problem.  To  incorporate  such  features  as  an  irregular  coast- 
line, irregular  bathymetry,  islands,  and  arbitrary  wind  stress  patterns, 
a  numerical  approach  becomes  necessary.  The  actual  numerical  approach  or 
algorithm  is  frequently  referred  to  as  a  numerical  model.  Models  may  be 
roughly  categorized  as  one-,  two-,  and  three-dimensional.  One-  and  two- 
dimensional   models  are  presented  in  this  chapter. 

The  one-dimensional   model  considers  motion  in  only  one  direction, 
e.g.   along  a  perpendicular  to  the  coastline.     It  is  presented  here  for  tv;o 
reasons:     its   simplicity  allows  physical    insight  into  the  problem  and  it 
provides  a  check  of  the  two-dimensional   model,  which  can  approach  a 
one-dimensional   model   as  a   limiting  condition. 

The  development  of  the  one-dimensional   numerical  model   begins  by 
referring  to  Equations  2-6. 


^  .  2co(sin<^)q^  =  -  ^  %  .  go  I  .  1  ir^^  -  x^^)  (2-6b) 
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^+^+1^=0  (2-6c) 

ax        dy        dt  ^         ' 

Considering  only  motion  in  the  x-direction  leads  to  Equations   (3-la)  and 
(3-1 b) 

A  A 

Using  the  quadratic  friction  law,  x,  =  — on2 —  >  Equations 

X 

(3-la)  and  (3-1 b)  can  be  written  as 

8q 

ir  ^  It  ="  t3-2b) 

These  are  the  equations  to  be  used  for  the  one-dimensional  calculations. 

The  finite  difference  integration  scheme  used  for  the  numerical 
calculations  is  similar  to  that  of  Reid  and  Bodine  [10]  and  Verma  and 
Dean  [28].  Figure  3-1  provides  a  schematic  diagram  of  a  typical  coastal 
region  showing  the  coordinate  system  and  defining  the  variables  in  use. 
The  vertical  coordinate,  z,  is  oriented  positively  upward  and  referenced 
to  mean  sea  level,  MSL;  the  horizontal  coordinate  x  is  oriented  toward 
shore  with  x  =  0  occurring  at  the  shelf  break.  In  the  computations  n> 
h,  T,  and  p  are  considered  uniform  over  the  grid  width  ax, and  centered 
at  the  midpoint  of  the  segment.  The  quantity,  q  ,  denotes  tne  flux  in 

A 

volume  per  unit  width  through  the  plane  representing  the  seaward  bound- 
ary of  a  grid  element. 
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For  clarity  the  variables  are  placed  into  an   indicia!    form  where 
q   (i+l,n)   indicates  the  volume  transport  per  unit  width  through  the 
plane  section  at  the  location  x  =   i   ax  and  at  time  nAt,  i.e.    flow  into 
the  grid  designated  by  index  i+1.     Similarly,  n(i,n+l/2)   indicates  the 
value  of  n  for  the  i —  grid  element.      It  is  assumed  to  be  centered  at 

» 

X  =  (i+l/2)Ax  and  to  be  evaluated  at  time  t  =  (n+l/2)At. 

Following  the  scheme  in  Appendix  A,  the  equations  are  placed  in 
finite  difference  form  as  Equations  (3-3a)  and  (3-3b) 

{[n(i-l,n+l/2)  -  n(i,n+l/2)]  g  +  ^  [p^(i-l  ,n+l/2)  -  p^(i  ,n+l/2)]} 

+  f^[T   (i-l,n+l/2)  +  T   (i,n+l/2)]]  (3-3a) 

x  x 

where 

The  continuity  Equation  (3-2b)  provides,  in  finite  difference  form, 
the  equation  for  the  free  surface, 

n(i,n+3/2)  =  n(i,n+l/2)  +  |^  [q^Ci.n+l)  -  q^(i+l,n+l)]         (3-3b) 

To  illustrate  the  calculation  procedure  assume  that  all  values  of 
q  are  available  up  to  and  including  t  =  nAt.  Beginning  at  the 
seaward  boundary,  calculate  q  's  at  time  (n+l)At  using  Equations  (3-3a). 

A 

Next,  calculate  n's  at  (n+3/2)At  using  Equation  (3-3b),  then  advance  the 
pressure  and  shear  stress  to  (n+3/2)At  and  calculate  q  's  at  time 
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(n+2)At,  and  so  forth.  In  effect,  averages  over  small  segments  are  being 
made  in  both  space  and  time  to  approximate  the  true  solution.  This  re- 
sults in  a  type  of  leapfrog  method,  since  q  's  at  time  (n+l)At  are  used 

A 

to  extend  the  solution  from  n(i,n+l/2)  to  ri(i.n+3/2)  and  similarly 
n  and  D  at  (n+3/2)  are  used  to  obtain  q  (i,n+2).  The  process  is  con- 

A 

tinued  until  a  desired  time  is  reached,  such  as  a  surge  maximum  at  the 
coast. 

The  seaward  boundary  condition,  the  barometric  tide,  is  applied 
to  grid  (l,n),  where  the  wind  setup  is  zero  and  h(l,n)  ^  «>.  This  con- 
dition is  written  as  n(l  ,n)  =  1 .15  [p  -  p  (l,n)]  where  n  is  in  feet  and 

oo    r] 

p^  is  the  undisturbed  pressure  in  inches  of  mercury  far  from  the  storm. 
The  landward  boundary  condition  of  no  flow  requires  q  (ic,n)  =  0.  The 
initial  conditions  used  will  require  all  q  's  =  0  and  the  grid  to  be  in 

A 

barometric  equilibrium.  The  only  remaining  data  required  for  the  calcu- 
lations are  the  pressure  and  shear  stress  distributions  as  functions  of 
time.  Figure  3-2  shows  the  movement  of  such  a  system  at  constant 
velocity  and  with  constant  form.  The  required  wind  and  pressure  data  are 
extrapolated  linearly  from  the  storm  model  described  in  Section  2-c. 

The  model  also  allows  for  the  "flooding"  or  "draining"  of  a  grid 
element.  Consider  the  element  (ic,n)  in  Figure  3-1,  if  ri(ic-l,n) 
<-h(ic)  then  the  system  is  as  shown,  if  n(ic-l,n)  >  -  h(ic)  then  the  grid 
element  (ic)  must  be  flooded  or  activated.  In  practice  when  activating 
a  grid  element  an  instability  can  occur  which  results  in  switching  be- 
tween the  activated  and  inactivated  states.  Also  remember  that  at  the 
moment  of  grid  activation  D  =  0  causing  trouble  with  the  calculations  at 
that  point.  To  avoid  these  problems,  the  following  steps  are  taken. 
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Once  activated  the  grid  element  is  constrained  to  remain  activated  fo.r  a 
specified  number  of  time  steps.  Five  time  steps  has  been  found  satis- 
factory in  the  present  model.  In  addition  to  the  above,  D  is  set  equal 
to  .01  ft.  in  the  just-flooded  grid  element  prior  to  the  calculations  for 
the  next  time  step  to  prevent  zero  divide  problems.  This  introduces  a 
small  error  into  the  system.  Since  flooding  occurs  infrequently  with 
respect  to  the  number  of  time  steps  this  error  is  considered  acceptable. 
To  calculate  a  water  velocity  u  from  q  an  average  depth  over  ad- 

A 

jacent  grids  is  taken  and  u(i,n)  is  approximated  as 

u(i.n)  =  2q^(i,n)  [D(i-1 ,n+l/2)  +  D(i ,n+l/2)]"'  (3-4) 

Since  u(i  ,n)  is  not  used  in  the  calculations,  the  approximations  in 
Equation  (3-4)  are  not  considered  significant.  A  more  accurate  value, 
centered  at  time  n,  which  must  be  calculated  one  time  step  after  the 
q  's  is  given  by  Equation  (3-5). 

A 

u(i,n)  =  4q^(i,n)  [D(i-1 ,n-l/2)  +  D(i-l,n+l/2)  +  D(i,n-l/2) 

+  D(i.n+l/2)]"'  (3-5) 

It  1s  assumed  that  the  hurricane  or  pressure  and  shear  stress  dis- 
tributions approach  normal  to  the  coast.  Numerical  stability  for  the 

1  /  O 

one-dimensional  model  requires  that  t—  <  (gD^^  )    where  D^^^  is  the 

AX  —   max  max 

maximum  total  depth  expected  anywhere  in  the  system.  Thus  given  a  bottom 
profile  and,  shear  stress  and  pressure,  as  a  function  of  time,  the  surge 
can  be  calculated  using  this  reasonably  simple  technique.  This  model  is 
particularly  useful  in  quantifying  effects  of  the  numerical  model  that  are 
not  inherent  in  the  equations  of  motion  but  result  as  an  effect  of  dis- 
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cretlzation.  The  model  is  also  easily  checked  and  thus  provides  a  basis 
for  comparison  with  other  models  which  are  not  easily  checked. 
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b.  Development  in  Two  Dimensions 

The  development  in  two  dimensions  is  analogous  to  that  in  one 
dimension  and  therefore  much  of  the  detail  will  be  omitted.  Recalling 

that  the  two-dimensional  model  considers  motion  in  the  x  and  y  directions, 

refer  to  Equation  (2-6),  and  employing  the  quadratic  bottom  friction  obtain 

3^  -  2.(sin,)q^  .  go  |a  =  .  ^  ^  .  -x  .  __^       (3-6a) 


!i...(3in,)q,.,0|a=-^!!..>.^ 


These  equations  are  now  in  forms  which  can  be  applied  and  where  no 
other  approximations  need  be  made.     Expressing  Equation   (3-6a)   in  finite 
difference  form  as  described  in  Appendix  B  leads  to  the  following  for  the 
x-equation. 

qx(i.J.n+l)  =   p-  [q^(i,j,n)  +  2w(sin4))q^  (i,j,n)At 


+  f^{T       (i-l.j,n+l/2)  +  T       (i,j.n+l/2)} 
^         x  X 

+  f— -  {D(i-l.j,n+l/2)  +  D(i.j,n+l/2)}   {  1  [p^(i-l  ,j,n+l/2) 

-P^(i  J.n+V2)]  +  g  [n(i-l  J  ,n+l/2)   -  ri(i  J  ,n+l/2)]}]  (3-7a) 
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where 


F^  =   1   +  ^  {   [q^(i,j,n]2  +  [qy(i,j,n)]2}''/2 


[D(i-l,j.n+l/2)  +  D(iJ,n+l/2]-2  (3-7b) 

Equations   (3- 7a),   (3-7b),    (3-8a),  and   (3-8b)   are  adapted  using  the 
form  given  by  Reid  and  Bodine  [10],  with  changes  in  notation. 
Following  the  same  procedure  for  the  y-equation  yields 

qy(i,j,n+l)   =  ■^[qy(i,j,n)   -  2oj(sin(())  q'^(i,j,n)At 
1^  {t     (i,j,n+l/2)  +  T      (i,j-l,n+l/2)} 

y  y 

+  ^  {D(i,j,n+l/2)  +  D(i,j-l,n+l/2)}     {i  [p^(i  ,j-l  ,n+l/2) 


P^(i,J,n+l/2)]  +  g[n(i,j-l,n+l/2)   -  n(i  ,j  ,n+l/2)]}]  (3-8a) 


where 


Fy  =  1   +  ^  {[qx(i,j.n)]2  +  q   (i,j,n)]2}  ^^^ 


[D(i,j,n+l/2)   +  D(i,j-l,n+l/2)]-2  (3-8b) 

Here,  ax  and  Ay  are  the  grid  spacings  as   shown   in   Figure  3-3.      Figure 
3-3  also  defines  the  variables   in  Equation   (3-7a)   and  the  sense  of 
operation  of  a  grid  element. 

The  values   for  the  variables,  n,  p   ,  h,  D,  t  are  assumed  to  be 

n 

centered  in  each  grid  as  shown  in  Figure  3-3.  The  values  are,  after 
that,  assumed  uniform  over  the  entire  grid  element  for  the  calculation 
at  that  time  step. 

Since  q^(i,j,n)  and  q  (i,j,n)  occur  at  different  points  it  becomes 
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necessary  to  obtain  spatially  averaged  values  q„(t,j,n)  and  q  (i,j,n) 

X  y 

for  use  in  the  Coriolis  term  and  the  friction  terms  F  ,  F  .  Thus, 

X   y 

q,,(i.j>n)  and  q'^(i,j,n)  are  considered  to  be  at  the  same  point  and 

qj^(i.j.n)  and  q"  (i,j,n)  are  considered  to  be  at  the  same  point.  This 

is    illustrated  more  clearly  in  Figures  3-4a  and  3-4b.  q"  Ci.j.n) 
is  defined  as 

qy(i,j,n)  =  |[qy(i,j,n)  +  qyCi-l,j,n)  +  q^(i,j+l,n)  +  q^(i-l  ,j+l  ,n)] 

and  q     is  expressed  by 

qx(i.J.n)   =|[q^(i,j,n)  +  q^(i+l.j,n)  +  q^O+l  ,j-l  ,n)  +  q^(i,j-l,n)] 


The  equation  of  continuity.  Equation  (3-6c),  becomes  in  finite 
difference  form, 

n(i,j,n+3/2)  -  n(i,J,n+l/2)  +  f^  [q^(i-l  ,j  ,n+l  )-q^(i  ,j  ,n+l )] 

+  fy  [qy(i.j-l,n+l)  -  q^(i,j,n+l)]  (3-9a) 

The  total  depth  D  is  defined  by 

D(i,j.O  =  h(i,j)  +  n(i,j,c)  (3-9b) 

where  5  is  an  arbitrary  time  increment. 

To  carry  out  an  actual    computation,  a  grid  network  is  set  up 
as  illustrated  in  Figure  3-5. 

The  coordinate  convention   is  the  same  as   Figure  2-1   with  x- 
positivG  toward  shore  and  y-positive  to  the  left  while  facing  the 
shoreline.     The  boundary  conditions  are  illustrated  in  Figure  3-5  and 
are  described  as   follows:     1)     the  fluid  flux  across  any  sea-land 
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FIGURE    3-5        SCHEMA     OF    GRID    SYSTEM     AND 

BOUNDARY     CONDITIONS     EMPLOYED     IN 
TWO  -  DIMENSIONAL     CALCULATIONS 
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boundary  is  equal  to  zero;  2)  the  boundary  condition  at  the  seaward 
boundary  is  given  as  n  equal  to  the  barometric  tide  plus  the  astro- 
nomical tide  at  that  point;  3)  at  the  lateral  boundaries  q  =  0; 
4)  also  at  the  lateral  boundaries  n  is  set  equal  to  the  barometric 
tide  plus  the  astronomical  tide.  Furthermore  it  will  be  assumed  that 
the  lateral  boundaries  are  far  enough  away  from  the  region  of  interest 
that  the  lateral  boundary  conditions  have  little  effect.  It  should  be 
reemphasized  that  the  lateral  conditions  are  those  conditions  arrived 
at  empirically,  which  allow  calculations  without  instabilities  or  other 
problems  with  the  model.  Initially  the  system  is  considered  quiescent 
and  in  barometric  equilibrium.  This  condition  is  of  little  effect 
providing  the  storm  begins  sufficiently  far  offshore.  In  this  model 
the  bathymetry  and  shoreline  are  arbitrary  and  the  presence  of  islands 
is  allowed. 

Currents  in  the  coastal  waters  are  of  interest  for  several 
reasons,  among  them:  1)  to  aid  in  making  and  evaluating  future  current 
and  surge  calculations,  2)  at  the  coastline  hurricane-generated 
currents  are  of  interest  as  a  cause  of  beach  erosion,  3)  hurricane- 
generated  currents  may  be  instrumental  in  producing  extreme  forces  on 
ocean  structures,  and  4)  currents  are  instrumental  in  dispersion  and 
pollution  studies.  Thus,  it  is  of  further  interest  to  describe  the 
scheme  used  to  calculate  the  water  velocity  at  a  specified  grid 
element.  As  in  the  one-dimensional  case  the  total  depth  at  time  n-1/2 
is  used  to  calculate  u(i,j,n)  and  v(i,j,n).  In  order  to  maintain  a 
consistent  approach,  an  average  in  both  the  x  and  y  directions  is 
taken.  This  leads  to  the  following  value  of  u(i,j,n)  and  for  v(i,j,n) 
as  described  in  Figure  3-6. 
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u(i.j.n)  =  8q^(i,j,n)   [2D(i-l  ,j  ,n)+2D(i  ,j  ,n)+D(i  ,j+l  ,n)+D(i-l  .j-1  ,n) 
+D(i.j-l.n)+D(i-l,j+l,n)]"> 

v(i,j,n)   =  8q^(i,j,n)    [2D( i  ,j-l  ,n)+2D(i J  ,n)+D(i-l ,j ,n)+D(i-l ,j-l ,n) 
+D(i+l,j,n)+D(i+1,j-l,n)]"i 

Adjacent  to  the  land-sea  boundary  it  is  necessary  to  modify  this 
procedure  because  of  the  no  flow  boundary  condition. 
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c.  Effects  of  Discretization  and  Checks  of  the  Numerical  Model 

In  any  model  of  this  type  it  is  desirable  to  obtain  some 
estimate  of  the  effects  of  the  grid  size,  time  increments,  and  lateral 
boundary  conditions  on  the  model  output.  A  check  of  the  programming 
is  also  necessary.  One  way  to  accomplish  this  is  with  a  series  of  two 
or  more  programs  which  although  of  different  concept  or  design,  converge 
to  identical  solutions  for  some  limiting  case.  This  does  not  neces- 
sarily prove  all  are  correct  but  provides  confidence  in  the  results. 

A  one-dimensional  time  dependent  program  should  give  results 
similar  to  those  of  a  two-dimensional  time  dependent  program  if  the 
radius  scale  of  the  hurricane  is  large  enough  with  respect  to  some 
length  scale  characteristic  of  the  system.  A  comparison  of  this  type 
has  been  carried  out  for  an  extremely  large  storm,  with  the  results 
and  storm  characteristics  as  shown  in  Figure  3-7.  One  would,  in 
general,  expect  the  one-dimensional  program  to  over  predict  water 
levels  far  from  the  coast  since  a  cross-section  of  the  storm  is  used 
and  this  effect  extends  laterally  to  infinity,  but  in  the  two-dimen- 
sional case  the  pressures  and  stresses  die  off  with  distance  from  the 
center  line.  Near  the  coast  the  Coriolis  effect  and  the  two-dimensional 
effects  dominate.  The  one-dimensional  plot  shown  in  Figure  3-7a  repre- 
sents a  cross  section  of  maximum  surge  heights.  The  bottom  was  a 
linear  slope  approximated  by  a  finite  difference  form. 
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As  a  further  comparison,  a  time  independent  or  steady  state 
model  was  developed  as  a  further  check  of  the  computer  modeling.  It 
is  one  dimensional  and  time  independent  in  that  irrr-    is  equal  to  zero, 
q  must  then  be  equal  to  zero  and  the  bottom  friction  term  is 
eliminated.  The  numerical  technique  consists  of  a  Runge-Kutta 
starter  and  Adam's  Method  predictor  [29].  The  stepwise  bottom  profile, 
with  slope  .0005  as  before,  was  included  and  flooding  allowed,  also 
as  before.  The  storm  characteristics  for  the  one-dimensional  time 
independent  case  are  shown  in  Figure  3-7b,  c,  with  the  storm  center 
placed  10  nautical  miles  offshore.  Note  that  although  the  storm  is 
stationary  it  has  the  characteristics  of  a  storm  traveling  with 
V  =  10  knots.  The  numerical  results  and  the  steplike  bottom  profile 
are  illustrated  by  the  dotted  line  in  Figure  3-8. 

In  order  to  compare  these  results  with  those  of  the  one- 
dimensional  time-dependent  routine  the  time-dependent  routine  was 
run  with  identical  storm  characteristics  and  then  stopped  when  the 
storm  center  was  10  nautical  miles  offshore.  In  this  case  the 
bottom  friction,  f,  was  set  equal  to  zero  and  as  a  result  the  solution 
oscillated.  The  shaded  boxes  in  Figure  3-8  represent  the  range  of 
oscillations.  The  stepwise  approximation  to  the  uniform  slope, 
.0005,  was  again  used  in  this  case. 

How  well  does  the  10  nautical  mile  step-like  bottom  profile 
model  a  true  linear  slope?  This  question  can  be  answered  as  follows. 
The  time  independent  model  was  run  with  all  variables  as  before,  with 
the  single  exception  of  a  uniformly  sloping  bottom.  These  results  are 
plotted  as  the  solid  line  in  Figure  3-8.  In  deep  water  the  results  are 
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nearly  identical.     In  the  shallow  region  the  differences  become 
noticeable.     At  the  risk  of  over  generalizing,  this  stepwise  grid 
network  is  permissible  for  values  of  D  =  h  +  n  greater  than  20  feet; 
for  values  less  than  20  feet  the  errors  can  become  significant. 

In  the  two-dimensional   problem  another  "error-dimension" 
exists,  that  of  the  effect  of  placement  of  lateral   boundaries.     This 
is,  in  addition  to,  the  effects  of  discretization.     The  lateral   bound- 
aries of  the  model   are  an  artifice     and  if  they  are  not  a  sufficient 
distance  from  the  zone  of  high  wind  stress     then  their  proximity  will 
influence  the  calculated  storm  tides  and  velocities.     In  order  to 
evaluate  this  effect  some  basis  fo>"  comparison  must  be  available.     A 
two-dimensional   analytical   treatment  of  storm  tides  on  continental 
shelves  is  available  from  Dean  and  Pearce  [21]. 

Consider  the  situation   illustrated  in   Figure  3-9  in  which  a 
"strip"  wind  field  is   uniform  within  a  total   width,  w,  and  zero  out- 
side of  this  width.     The  Coriolis  effect  is  omitted  and  the  bottom 
shear  stress  is  expressed  as  proportional   to  the  transport  components 
(i.e.   linearized).     Only  shear  stresses  in  the  horizontal    plane  are 
considered,  i.e.   adjacent  water  columns  are  considered  to  be  laterally 
uncoupled  by  stress  terms;   it  may  be  anticipated  that  this  will   result 
in  a  discontinuity  in  q   ,  the  transport  component  in  the  x-direction. 

A 

The  windstress  in  the  region  |y|  <  w/2  drives  the  water  onshore  and 
results  in  a  storm  tide  field  which  is  a  maximum  at  (x'  =  x/2.  =  1.0, 
y'  =  y/ii  =  0).  The  q  field  in  this  region  is  a  maximum  at  the 
shelf  break  and  decreases  to  zero  at  the  coastline.  It  is  clear,  on  an 
intuitive  basis,  that  the  narrower  the  wind  field  relative  to  the 
shelf  width,  the  smaller  the  storm  tide  will  be.  For  a  very 
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broad  wind  field,  the  storm  tide  at  the  shelf  should  approach  the  value 
corresponding  to  a  one-dimensional  condition.  For  wind  fields  which 
are  laterally  limited,  there  will  be  a  gradient  in  the  storm  tide 
field  such  that  the  storm  tide  approaches  zero  as  |y'|  -»-  <=.  The 
surface  gradient  and  wind  stress  fields  will  cause  an  onshore  transport 
in  the  wind  stress  region,  an  offshore  transport  outside  of  this 
region  and  a  y  component  of  transport  ("  longshore  transport")  which 
is  directed  outward  from  the  region  of  onshore  wind  stress.  Response 
in  terms  of  dimensionless  n,  q  and  q  are  presented  in  Figures 
3-lOa,  b,  c. 

If  the  (artificial)  lateral  boundaries  in  the  numerical  model  are 
not  a  sufficient  distance  from  the  zone  of  high  wind  stress,  the 
proximity  of  the  model  boundaries  will  influence  the  calculated  storm 
tides  and  velocities.  In  order  to  evaluate  this  proximity  effect, 
numerical  model  computations  were  carried  out  for  the  shelf  geometry 
shown  in  Figure  3-11.  The  friction  and  total  water  depth  terms  in 
the  numerical  model  were  linearized  in  order  that  the  differences  in 
quantities  calculated  by  the  numerical  model  and  the  analytical 
solution  should  be  indicative  only  of  artificial  effects  introduced 
by  the  numerical  model.  The  results  of  these  calculations  are  de- 
scribed in  the  following  paragraphs. 

The  first  set  of  calculations  was  carried  out  for  a  ratio  of 
shelf  width  to  one-half  the  wind  field  breadth  (a  5  jl/(w/2))  of  1.83 
and  a  grid  size,  Ax,  which  is  variable.  The  lateral  extent  of  the 
shelf,  incorporated  into  the  model  is  denoted  w*  and  the  ratio  of 
wind  field  width  to  shelf  extent,  a*(i.e.  a*=w/w*).  The  parameter 


81 


t 


Nofes: 
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(3)  Results     Shown      Apply 
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a*  was  varied  over  the  range  0.24  <  a*  <  0.73.  The  results  [30] 

are  presented  in  Figure  3-12,  where  the  ordinate  represents  the 

maximum  storm  tide  ratio  as  determined  by  the  numerical  and  analytical 

solutions.  Perfect  agreement  is  indicated  by  a  value  of  1.0. 

From  this  figure,  if  a*  <  0.5,  the  resulting  errors  in  n   will  be 
^  ^         max 

less  than  10%.  As  indicated  in  Figure  3-12,  various  ratios  of 
shelf  length,  z,   to  grid  size  v/ere  employed,  however  there  does  not 
appear  to  be  any  significant  differences  in  the  results  for  2,/Ax 
values  as  low  as  3-1/2.  It  is  noted  that  these  calculations  were 
conducted  for  a  =  1.83  and  that  for  a  yery   wide  shelf,  the  parameter 
a**  =   x,/w*;  may  also  be  of  importance. 
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EFFECT     OF     LATERAL     BOUNDARY 
PROXIMITY      AND     RELATIVE     GRID     SIZE 
ON      MODEL     RESPONSE 


CHAPTER  4 
RESULTS  OF  NUMERICAL  CALCULATIONS 
a.  Introduction 

The  surges  and  currents  produced  by  Hurricane  Camille  were  chosen 
for  analysis  in  this  study.  The  severity  of  the  storm  and  the  extensive 
damage  resulting  from  it  crystallized  concern  about  storm  surge  calcu- 
lations. The  only  available  open-coast  current  data  taken  during  a  hur- 
ricane were  by  Murray  [16,17]  during  Hurricane  Camille.  The  surge  data 
available  for  hurricane  Camille,  consisting  primarily  of  high  water 
marks  [31,32]  while  by  no  means  extensive,  are  significant  when  compared 
to  data  from  other  storms.  Affected  by  the  storm,  is  the  section  of  the 
U.  S.  Gulf  Coast  ranging  from  Barataria  Bay,  Louisiana  to  Panama  City, 
Florida.  This  region  is  shown  in  Figure  4-1,  where  the  water  depth  con- 
tours are  in  fathoms.  For  the  purposes  of  the  numerical  calculation  the 
region  is  divided  into  a  grid  element  network  as  in  Figure  3-5.  Two 
grid  dimensions  were  used,  one  with  a  16.16  nautical  mile  grid  width, 
the  other  5.99  nautical  miles.  Square  grid  elements  are  used  in  both 
cases.  The  large  (16.16  n.m.)  grid  array  is  shown  in  Figure  4-2a 
superimposed  on  a  map  of  the  region,  the  grid  indexes  i,  j  are  given 
along  the  right  hand  side  and  bottom  of  the  figure,  respectively.  The 
depths  are  shown  in  Figure  4-2b  where  dry  land  is  indicated  by  a  minus 
number  and  the  location  of  the  decimal  point  represents  the  center  of  a 
grid  element.  The  depth  shown  indicates  an  average  depth  for  the  grid 
element  and  is  assumed  to  be  located  at  the  center  of  the  grid  element. 
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Figure  4-2c  shows  the  small  grid  model  and  Figure  4-2d  the  depths  used 
in  the  small  grid  model.  The  grid  size  and  arrangements  were  chosen 
heuristically  to  approximate  the  land  mass  as  closely  as  possible.  The 
storm  is  idealized  by  the  storm  model  discussed  in  Section  2-c.  The 
storm  parameters  are  assumed  constant  during  the  traverse  over  the  shelf. 
Therefore,  a  grid  network  of  storm  parameters  is  set  up  and  superimposed 
on  the  existing  grid  network.  This  network  of  storm  parameters  is  then 
moved  ashore  at  constant  velocity  with  the  values  of  surface  pressure 
and  shear  stress  interpolated  to  the  centers  of  the  stationary  grid 
elements.  Hurricane  Camille  is  assumed  to  travel  at  13  knots  and  to  be 
in  the  vicinity  of  Pass  Christian,  Mississippi  at  1130  hours  CDT  August 
17,  1969.  Also,  for  the  purpose  of  computations  it  is  assumed  that  the 
storm  motion  is  parallel  to  the  grid  lines  crossing  the  coast.  The 
idealized  and  actual  storm  tracks  are  shown  in  Figure  4-2a. 
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b .  Numerical  Calculations  and  Comparison 
with  Available  Data 

Figures  4-3a-e  depict  the  calculated  surge  heights  at  different 
times,  using  the  large  grid  model  and  storm  characteristics  for  Hurri- 
cane Camille.  Although  Figures  4-3a-e  provide  output  from  the  model  in 
an  interesting  form,  unfortunately  two-dimensional  arrays  of  measured 
surge  height  elevation  are  not  available  for  comparison.  In  order  to 
gain  some  idea  of  the  effectiveness  of  the  model,  comparisons  must  be 
made  with  measured  data.  This  section  will  emphasize  these  comparisons 
and  interpretations  to  be  made  from  them. 

Figure  4-4  depicts  the  measured  maximum  surge  that  occurred  along 
the  Gulf  Coast  during  Hurricane  Camille.  These  data  are  from  high  water 
marks  [31]  and  it  should  be  noted  that  a  maximum  high  water  mark  of  24.6' 
was  determined  at  Pass  Christian,  Mississippi  from  a  "debris  line". 
Superimposed  on  Figure  4-4  are  the  numerical  calculations  using  the 
previously  stated  assumptions.  The  values  shown  represent  maximum 
surges  at  all  locations.  A  bottom  friction  coefficient  of  f  =  .005  was 
used  for  the  calculations  shown  in  Figure  4-4.  It  can  be  seen  in 
Figure  4-5  however,  that  the  calculated  surge  is  not  sensitive  to  f. 
The  model  used  to  obtain  Figure  4-5  was  an  earlier  version  than  the 
present  model,  and  did  not  include  grid  elements  (15,9),  (16,9),  (15,7), 
(16,7),  and  (15,6).  The  sensitivity  of  the  model  is,  however,  clearly 
illustrated. 

In  general  the  locations  of  calculated  surges,  centered  at  the 
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grid  elements,  will  not  coincide  with  actual  coastline  locations  or 
other  locations  of  interest.  In  order  to  establish  surge  elevations  at 
locations  other  than  at  those  calculated,  two  techniques  were  used.  To 
obtain  the  surge  at  a  location  not  on  the  coast,  such  as  an  island,  a 
direct  correspondence  between  model  and  prototype  was  assumed,  as  shov/n 
in  Figures  4-2a  and  4-2c.  To  obtain  the  surge  at  the  location  in 
question,  the  values  at  the  surrounding  grid  centers  are  used  to  find  a 
value  at  the  desired  point  by  linear  interpolation. 

Since  the  effects  of  the  boundary  may  be  important,  a  different 
scheme  is  used  to  correlate  between  points  on  the  model  coastline  and 
the  prototype  coastline.  If  the  grid  boundaries  and  the  coastline  coin- 
cide there  is  no  problem.  When  the  coastline  and  grid  boundary  do  not 
exactly  coincide  then  the  point  on  the  model  coastline  nearest  the  point 
on  the  prototype  coastline  is  chosen  as  the  point  at  which  to  obtain  the 
surges  and  currents.  This  technique,  could  lead  to  inconsistencies;  how- 
ever, all  calculations  for  a  particular  comparison  are  done  in  the  same 
manner.  In  general  since  the  available  data  are  from  high  water  marks 
it  is  necessary  to  extrapolate  to  the  coast  as  shown  in  Figure  4-6. 


7?(i,  j) 


JlillL 


i 


mill 


,17      =•7(1  +  1,  j)  +  17(14-1,  j)--n(i,  j) 

COAST  — = ■ 


FIGURE  4-6  SCHEME  USED  TO  EXTIV.POLATE  SURGE  HEIGHTS  TO  COAST 
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Because  the  storm  surge  is  not  linear  with  distance,  linear  extrapo- 
lation produces  an  error.  Since  the  model  is  discrete  with  a  varying 
forcing  function,  it  is  unreasonable  to  extrapolate  over  many  grid  ele- 
ments; for  this  reason,  and  also  because  of  its  simplicity,  linear  extra- 
polation is  used.  Figure  4-7  compares  the  calculated  tidal  heights  at 
the  mouth  of  Biloxi  Bay,  Mississippi  for  both  the  extrapolated  and  not 
extrapolated  cases.  Grid  elements  (8,  13)  and  (9,  13)  are  used  to 
extrapolate  to  the  coast. 

The  grid  size  of  the  model  is  too  large  to  appropriately  handle 
small  bays  and  estuaries.  The  barrier  islands  in  the  region,  such  as 
Dauphine  Island,  Alabama,  have  a  characteristic  long  thin  shape  which 
presents  a  problem  in  modeling  them.  Their  area  is  small  compared  to  a 
grid  element,  making  it  unreasonable  to  use  whole  grid  elements  for  the 
island,  yet  they  may  provide  an  effective  flow  barrier  over  the  width  of 
one  or  more  grid  elements.  Bearing  in  mind  that  these  islands  were  in 
general  flooded  during  the  hurricane  passage  the  problem  was  approached 
as  follows:  For  the  large  grid  model  the  islands  were  ignored.  For  the 
small  grid  model  the  flows  at  appropriate  grid  boundaries  were  constrain- 
ed to  be  zero  and  for  comparison  the  small  grid  model  was  also  run  with 
no  flow  restrictions.  In  all  of  these  cases  a  bottom  friction  coeffi- 
cient of  f  =  .005  was  used. 

Figure  4-8a  depicts  the  measured  water  level  at  the  Dauphine 
Island,  Alabama  Marine  Laboratory  [31]  compared  to  the  surge  calculated 
by  the  large  grid  model  and  interpolated  to  the  appropriate  location. 
There  are  very  little  tide  data  available  and  the  Dauphine  Island  data 
are  unique  in  being  close  to  the  open  Gulf.  Most  tide  data  are  taken  in 
rivers  and  estuaries  where  the  tidal  conditions  are  different  from  those 
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in  the  Gulf  of  Mexico.  Also  superimposed  cri  these  data  is  the  predicted 
astronomical  tide  [33]. 

Figure  4-8b  depicts  the  measured  tide  at  Dauphin  Island  compared 
to  results  of  the  small  grid  model  with  and  without  the  flow  restriction 
representing  the  island.  Note  the  difference  between  the  Gulf  side  and 
the  Mississippi  Sound  side  of  the  island  for  the  calculations  with  flow 
restriction.  Note  also  the  increased  phase  difference.  A  weighted 
average  of  the  elements  (23,14) ,  (23,13) ,  (22,14),  and  (22,13)  of  the 
small  grid  model  was  used  to  locate  the  calculations  at  Dauphine  Island. 
The  effect  of  Mobile  Bay  which  is  included  in  the  small  grid  model  may 
produce  the  discrepancies  noted  above.  The  assumption  of  constant 
storm  characteristics  may  also  be  responsible  for  the  phase  difference 
of  measured  and  predicted  surge  height. 

Average  water  currents  were  also  calculated.  The  current  measure- 
ments by  Murray  [16,17]  are  plotted  along  with  direction  and  wind 
measurements  in  Figure  4-9a.  Figure  4-9b  shows  the  coastal  topography 
in  the  immediate  vicinity  of  the  current  measurements.  Notice  that  the 
beach  is  oriented  East-West  and  that  at  the  time  of  maximum  currents  the 
flow  direction  is  offshore.  The  meter  was  installed  1.4  meters  above 
the  bottom,  in  6,3  meters  of  water  and  about  360  meters  offshore. 

Figure  4-10  depicts  the  observed  current  magnitude  from  Reference 
[16]  as  compared  to  four  cases  of  f  and  X  for  the  large  grid  model  and 
one  case  for  the  small  grid  model.  The  most  important  point  here  is  the 
sensitivity  of  the  model  to  changes  of  both  bottom  friction  f  and  the 
starting  point  of  the  storm  X^,  in  nautical  miles  from  storm  landfall. 

Recent  longshore  current  measurements  [34]  were  made  approximately 
fifteen  hours  after  the  passage  of  Hurricane  Agnes  (1972).  These 
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results,  presented  in  Figure  4-11,  show  that  longshore  currents  of  2-3 
ft/sec.  are  entirely  reasonable.  The  beach  slope  of  m  =  .02  is  the  same 
as  for  the  condition  at  Eglin  AFB  and  the  wave  period  of  8  sec.  slightly 
less  than  the  10  sec.  maximum  observed  by  Sonu  [35].  Of  interest  is  the 
magnitude  of  the  calculated  winds  for  Camille  and  a  comparison  with 
measured  data  in  Figure  4-12. 

In  examining  Murray's  data,  one  is  tempted  to  postulate  the  exist- 
ence or  even  predominance  of  longshore  currents,  however  it  is  not 
reasonable  to  explain  the  current  normal  to  shore  as  a  longshore  current. 
The  numerical  model  predicts  essentially  zero  normal  current  that  close 
to  shore  since  the  normal  flux  through  the  boundary  is  constrained  to 
zero.  Since  the  wind  direction  is  approximately  south  preceding  the  off- 
shore current  a  quasi -steady-state  circulation  cell  with  onshore  flow  at 
the  surface  and  offshore  flow  at  the  bottom  could  explain  the  observed 
phenomena.  The  important  point  is  that  the  assumptions  made  in  the  pres- 
ent numerical  model  preclude  the  representation  of  a  flow  variation  over 
depth.  An  obvious  conclusion  is  that  more  extensive  data  are  needed  to 
better  explain  the  processes  involved  in  nearshore  circulation. 

Damage  to  low-lying  coastal  areas  is  caused  primarily  by  flooding 
and  the  ensuing  wave  action  and  beach  erosion.  To  a  "user"  of  a  storm 
Surge  prediction  scheme,  frequently  the  only  result  of  importance  is 
whether  or  not  a  particular  area  will  be  flooded.  For  the  case  of  a  storm 
approaching  the  coast,  the  latest  available  storm  parameters  would  be 
used  in  predicting  flooding  to  warn  local  residents.  Design  of  various 
coastal  structures  relies  on  surge  and  flooding  calculated  for  a  statis- 
tically appropriate  storm,  frequently  the  Standard  Project  Hurricane,  SPH. 
A  detailed  study  of  flooding  is  not  within  the  scope  of  this  project. 
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LONGSHORE    CURRENT    AND    BEACH    PROFILE 
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However,  a  first  approximation  of  the  flooding  and  the  comparison  of 
measured  and  predicted  flooding  provides  a  suitable  conclusion  to  this 
stage  of  the  study. 

The  flooded  region  is  predicted  using  the  simple  assumption  that 
the  water  level  inland  near  the  coast  will  be  at  the  predicted  surge 
elevation.  Thus  the  water  is  assumed  to  extend  inland  horizontally  from 
the  elevation  predicted  at  the  coast.  For  example,  if  a  surge  10  feet 
above  MSLv/ere  predicted  for  a  region,  then  the  predicted  flooded  area 
for  that  region  would  extend  from  the  normal  coastline  to  the  10  foot 
Cabove  MSL)  contours  on  a  topographical  map.  The  predicted  surge  in 
this  case  is  that  from  the  large  grid  model  as  shown  in  Figure  4-4,  with 
the  surge  elevation  extrapolated  to  the  coastline.  A  typical  comparison 
of  these  predictions  and  the  actual  flooding  as  documented  by  the  U.  S. 
Army  Corps  of  Engineers  [31]  is  given  in  Figure  4-13. 

If  the  results  of  the  type  in  Figure  4-13  are  to  be  evaluated,  some 
quantitative  measure  of  the  "goodness"  of  the  fit  to  the  data  is  neces- 
sary. Because  of  the  convoluted  nature  of  the  flooded  region  comparisons 
based  on  a  flooded  length  would,  in  general,  have  little  significance. 
For  this  reason  it  was  decided  to  compare  actual  and  predicted  flooded 
areas  per  unit  length  of  coastline.  In  practice  this  consisted  of  taking 
an  approximately  3  mile  wide  strip  running  north  from  the  coast,  and,  for 
that  strip,  calculating  the  measured  and  predicted  flooded  areas  as  de- 
termined using  a  planimeter.  As  a  consistent  measure  these  areas  were 
divided  by  the  width  of  the  strip  to  obtain  area  per  unit  length.  When- 
ever possible  the  areas  were  taken  to  the  northern  extent  of  the  flooding. 
In  several  cases  this  was  impossible  and  an  arbitrary  boundary  was  drawn 
about  8  statute  miles  from  the  coast. 
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The  plot  in  Figure  4-14  runs  from  Edgev/ater  Park,  Mississippi 
eastward  to  Bayou  La  Batre,  Mississippi.  In  general  this  scheme  over- 
predicts  the  flooded  area.  There  are  two  apparent  reasons  for  this: 
1)  the  model  predicts  surges  slightly  too  large  in  that  region,  2)  fric- 
tion and  inertia  prevent  the  surge  from  reaching  points  far  inland  even 
though  they  are  at  a  low  elevation. 
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CHAPTER  5 
FIELD  DATA  AND  INSTRUMENTATION 

The  measurement  of  stom  surges  and  currents  is  difficult 
because  of  the  extreme  weather  conditions  prevailing  at  the  time 
observations  are  to  be  made.  There  exist  two  basic  methods  for 
approaching  this  problem,  excluding  that  of  the  observer  being  present 
during  the  storm.  One  may  have  a  permanent  installation  or  install- 
ations and  wait  for  a  significant  storm,  or  one  may  "run  before  the 
storm"  and  position  equipment  to  record  the  necessary  phenomena. 
The  latter  appears  the  best  choice  but  the  logistics  of  implementation 
and  recovery  are  great,  and  seas  may  be  quite  rough  preceding  the 
storm  due  to  the  presence  of  forerunners.  Sonu  [35]  taking  obser- 
vations at  Santa  Rosa  Island,  Florida  observed  a  period  of  wave 
action  dominated  by  forerunners  originating  in  Hurricane  Camille.  This 
period  lasted  for  28  hours  from  1200  hours  Central  Daylight  Time,  (CDT) 
August  16  to  1600  hours  CDT  August  17  and  ended  with  the  arrival  of 
gale-force  winds  and  locally  generated  windwaves.  The  period  of 
forerunners  was  characterized  at  this  location  by  a  maximum  breaking 
wave  height  of  approximately  3  meters.  At  the  closest  point  of 
approach  Hurricane  Camille  was  200  miles  from  the  point  of  observation. 

The  considerations  for  the  predictions  of  hurricane  movement  are 
not  within  the  scope  of  this  report.  The  predictions  are,  however. 
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an  important  consideration  for  successful  completion  of  any  field 
measurement  program.  To  acquire  insight  into  the  problems  involved 
in  gathering  field  data,  observe  the  hurricane  advisories  issued  by 
the  National  Weather  Service  as  presented  in  Figure  5-1  and  also 
observe  that  storm  landfall  near  Bay  St.  Louis,  Mississippi  between 
2300  and  2400  hours,  CDT,  on  August  17,  1969.  Using  this  information, 
a  Gedanken*  field  measurement  can  be  carried  out.  At  1900  hours  CDT 
the  storm  was  over  the  western  tip  of  Cuba,  traveling  to  the  north- 
west. At  that  time  the  field  crew  readied  all  instrumentation  and 
equipment.  The  point  of  landfall  is  unknovjn,  assuming  the  crew  to 
be  in  New  Orleans,  they  remain  since  the  probability  of  intercepting 
the  storm  does  not  increase  by  moving  at  this  time.  At  2300  hours 
CDT  August  15  the  National  Weather  Service  advised  that  a  hurricane 
watch  would  probably  be  established  for  the  northeast  Gulf  region 
prior  to  1100  hours  CDT  August  16.  At  this  time  the  crew  begins 
traveling  east  and  at  0800  hours  CDT  receive  the  bulletin  from  the 
National  Weather  Service  extending  a  hurricane  watch  from  Biloxi, 
Mississippi  to  St.  Marks,  Florida  (Figure  4-1).  It  is  decided 
to  install  the  instrumentation  at  the  midpoint  of  the  0800  watch 
area  near  Ft.  Walton  Beach,  Florida.  Three  hours  later  at  1100  hours 
Advisory  10  is  issued  posting  hurricane  warnings  from  Ft.  Walton  to 
St.  Marks.  By  this  time  the  instrumentation  has  been  installed  and 
the  further  success  is  left  to  chance.  The  hurricane  actually  went 
ashore  at  Bay  St.  Louis,  Mississippi  and  the  field  crew  missed  by 


*Gedanken  -  To  carry  out  by  proposing  a  hypothesis  in  thought  only. 
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Advisory  Number  5,  1100  CDT,  Friday  August  15. 

Cuba  warned  of  tides  up  to  8  ft.  and  5  to  10  inches  of  rain. 

Advisory  Number  6,  1700  CDT,  Friday  August  15. 

Gale  warnings  issued  for  Marquesas  Keys  and  Dry  Tortugas. 

Advi<=ory  Number  7,  2300  CDT,  Friday  August  15. 
Stated  Camille  a  dangerous  hurricane  and  poses  a  great  threat  to 
the  U.S.  mainland.  Advised  that  a  hurricane  watch  would  probably 
be  issued  for  a  portion  of  the  coastal  area  of  the  northeast 
Gulf  prior  to  1100  CDT,  Saturday  August  16. 

Advisory  Number  8,  0500  CDT,  Saturday  August  16. 
Small  craft  along  the  northwest  Florida  coast  and  as  far  west 
as  Mobile  urged  to  seek  safe  harbor.  Repeated  that  a  hurricane 
watch  v/ould  be  issued  later  in  morning. 

Advisory  Number  9,  0800  CDT,  Saturday  August  16. 
Hurricane  watch  issued  extending  from  Biloxi,  Mississippi,  to 
St.  Marks,  Florida.  Advised  that  a  hurricane  warning  would  be 
issued  at  1100  hours. 

Advisory  Number  10,  1100  CDT,  Saturday  August  16. 

Hurricane  warning  issued  extending  from  Ft.  Walton  to  St.  Marks. 

Gale  warnings  elsewhere  from  Pensacola  to  Cedar  Key. 

Advisory  Number  11 ,  1700  CDT,  Saturday  August  16. 

Highest  winds  raised  from  115  MPH  to  150  MPH.  5  to  12  feet 

tides  forecast. 

Advisory  Number  12,  2300  CDT,  Saturday  August  16. 
Tides  of  15  feet  forecast. 

Advisory  Number  13,  0500  CDT,  Sunday  August  17. 

Hurricane  warnings  extended  westward  to  Biloxi.  Hurricane 

watch  extended  westward  to  New  Orleans  and  Grand  Isle. 

Bulletin,  0500  CDT,  Sunday  August  17. 
Advised  that  if  present  movement  toward  the  mouth  of  the 
Mississippi  continued  for  a  few  more  hours,  hurricane  warnings 
would  be  extended  westward  to  New  Orleans  and  Grand  Isle. 

Advisory  Number  14,  0900  CDT,  Sunday  August  17. 

Hurricane  warnings  extended  westward  to  New  Orleans  and  Grand  Isle, 

Advisory  Number  16,  1500  CDT,  Sunday  August  17. 

Highest  winds  raised  to  190  MPH  and  tides  15  to  20  feet. 

Predicted  landfall  near  Gulfport,  Mississippi. 

FIGURE  5-1  SUMMARY  OF  CRITICAL  WARNINGS 
(FROM  REFERENCE  [36]) 
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200  miles.  They  did,  however,  recover  some  current  data,  which 

turned  out  to  be  the  only  data  in  existence  for  wind-generated  currents 

during  a  hurricane. 

It  should  be  noted  that  this  is  in  no  way  related  to  the  data 
of  Sonu  [35]  and  Murray  [16];  Murray's  instrumentation  was  in  place 
for  another  study  and  fortunately  current,  wind  data,  and  wave  height 
observations  were  obtained. 

This  hypothetical  deployment  of  equipment  in  anticipation  of  a 
hurricane  points  out  that  this  type  of  operation  involves  a  certain 
amount  of  risk,  since,  by  observation  the  storm  forerunners  were 
beginning  to  arrive  during  the  installation  procedure  (see  Figure 
5-2).  A  problem  postponing  the  operation  a  few  hours  may  postpone 
it  completely.  Thus,  matters  such  as  method  of  equipment  deployment 
are  very  important  and  in  fact  critical  since  the  operating 
environmental  conditions  will  be  extreme.  The  probability  of  loss 
of  equipment  is  impossible  to  determine  but  certainly  high.  Murray 
[16]  employed  three  current  meters,  two  in  the  nearshore  zone  which 
did  not  function  and  one  in  about  21  feet  of  water  seaward  of  the 
"outer  bar"  which  functioned  for  about  one  half  the  storm.  If  the 
storm  had  passed  directly  over  the  instrumentation  further  conse- 
quences could  have  been  expected.  A  system  that  is  100%  reliable 
is  desired  in  view  of  the  long  waiting  time  involved  and  in  view  of 
the  need  for  data  to  improve  surge  and  current  predictions. 

The  author  feels  that  the  technique  discussed  in  the  Gedanken 
measurement  is  not  reliable  enough.  A  system  that  is  deployable 
from  an  airplane  or  helicopter  would  be  very  desirable.  Also 
desirable  would  be  a  series  of  permanent  field  stations  located  in 
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FIGURE   5-2      HINDCAST    WAVE    FIELDS     FOR     HURRICANE 
CAMILLE       (AFTER    REFERENCE    [35]) 
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areas  of  high  hurricane  probability.  Harris  [37]  provides  in  Figure 
5-3  a  chart  of  absolute  number  of  destructive  storms  reaching  a  given 
location  in  a  54-year  period.  Using  this  and  other  data  including 
accessibility  of  the  area  in  question  it  would  be  desirable  to 
establish  a  network  of  five  permanent  underwater  surge  and  current 
monitoring  stations  around  the  State  of  Florida  in  the  configuration 
shown  in  Figure  5-4;  these  installations  would  be  activated  only 
during  the  hurricane  season.  This  network  in  conjunction  v/ith  a 
single  "run  before  the  storm"  installation,  when  practical,  should 
provide  the  necessary  redundancy. 

The  fixed  sites  should  have  a  recording  capability  of  30  days 
or  more.  This  would  then  require  three  or  four  servicings  of  the 
equipment  betvjeen  installation  at  the  beginning  of  the  hurricane 
season  and  recovery  at  the  end. 

The  development  of  inexpensive  but  effective  equipment  is  always 
desirable,  especially  where  a  high  probability  of  loss  is  involved. 
The  instrument  shown  in  Figure  5-5  and  Figure  5-6  is  a  recording 
underv/ater  tide  or  pressure  gauge  and  was  developed  from  items 
easily  available  commercially.  It  could  be  deployed  very  economically 
in  quantities  of  5-10  or  more.  At  the  present  time  the  response  time 
is  of  the  order  of  one  second.  For  the  prototype  situation  this  would 
be  increased  to  about  30  seconds  because  of  the  desirability  of 
suppressing  long  period  wave  action  present  during  a  storm.  A 
prototype  instrument  has  been  tested  at  Angel  fish  Creek  at  the 
northern  tip  of  Key  Largo  with  the  results  presented  in  Figure  5-7 
and  compared  to  a  nearby  conventional  tide  gauge  recorder. 
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FIGURE   5-3      NUMBER     OF     TIMES     DESTRUCTION 

WAS      CAUSED     BY    TROPICAL      STORMS. 
1901-1955       (FROM     REFERENCE     [37]) 
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Note: 


Instrumentation     Sites     Designated    by    ® 


FIGURE   5-4       PROPOSED     FIXED,     CONTINUOUSLY 

MAINTAINED    INSTRUMENTATION     SITES 
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FIGURE  5-6  RECORDING  UNDERWATER  TIDE  GAUGE 
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The  current  meter  presents  a  formidable  problem  since  the 
instrument  must  be  rugged  yet  sensitive  enough  to  obtain  a  mean  flow 
from  the  highly  turbulent  environment  existing  in  shallow  water  during 
storm  conditions.  The  current  meter  illustrated  in  Figures  5-8a  and 
5-8b  was  developed  in  conjunction  with  General  Oceanics,  Inc.*   It 
is  Inexpensive  but  proved  to  have  a  threshold  velocity  of  0.56 
feet/second,  high  enough  to  render  calibration  impossible  in  a  turbulent 
environment.  If  the  threshold  velocity  can  be  reduced  this  device 
would  prove  to  be  ideally  suited  for  the  measurement  of  hurricane 
driven  currents. 

Future  work  on  hurricane  surges  and  currents  is  especially 
needed  in  three  areas:  improved  representations  of  bottom  friction, 
considerably  more  field  data,  and  development  of  physical  modeling 
techniques.  It  is  evident  that  the  accurate  prediction  of  storm 
driven  currents  lies  in  a  better  understanding  of  the  dissipation 
mechanism  (i.e.  friction)  during  a  hurricane.  This  development  will 
rely  to  a  great  extent  on  availability  of  data  upon  which  a  better 
model  can  be  based.  Data  will  become  available  from  prototype 
situations  and  also  from  physical  models  of  which  there  are  none 
at  present.  There  are  many  problems  in  developing  a  physical  model. 
The  problem  of  scaling  is  difficult  since  the  length/depth  ratio  for 
a  typical  prototype  case  would  be  of  the  order  of  1,000-10,000.  Thus 
any  laboratory  model  would  probably  require  highly  distorted  length 
scales.  The  modeling  of  friction  resulting  from  turbulent  flows 
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FIGURE  5-8a     CURRENT  METER  DEPLOYED 
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FIGURE  5-8b     PHOTOGRAPHIC  RECORDING  MECHANISM 
OF  CURRENT  METER 
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and  pressure  and  shear  stress  distributions  will  all  require  a  sophis- 
ticated experimental  technique,  yet  would  produce  invaluable  data 
and  data  that  could  represent  a  considerable  equivalent  time  span  in 
the  prototype  situation. 


CHAPTER  6 
SUMMARY  AND  CONCLUSIONS 

The  response  of  coastal  or  shelf  waters  to  the  passage  of  an  orga- 
nized storm  system  is  studied  by  1)  simplified  analytical  techniques, 
2)  finite  difference  numerical  calculations  and  3)  comparison  of  the  nu- 
merical calculations  to  field  observations.  The  momentum  equations  for 
fluid  motion  and  the  continuity  equation  are  placed  in  turbulent  form  with 
the  instantaneous  velocities  being  replaced  by  time  mean  velocities  and 
with  the  laminar  stresses  replaced  by  Reynolds  stresses. 

These  equations  are  then  integrated  in  the  vertical,  eliminating  the 
z-equation,  with  the  result  that  all  velocities  represent  an  average  over 
the  depth.  Boundary  conditions  for  the  model  boundaries  are  presented, 
including  the  surface  and  bottom  shear  stresses.  Initially  the  model  is 
assumed  quiet  and  the  disturbance  to  begin  traveling  far  enough  from  the 
region  of  interest  to  make  the  effect  of  the  initial  conditions  negligible. 

A  simple  analytical  hurricane  model  which  calculates  the  gradient 
winds  and  takes  the  forward  velocity  of  the  hurricane  into  consideration 
is  used  to  predict  storm  winds  from  the  central  pressure  index  Ap,  the 
radius  to  maximum  winds  R,  and  the  forward  velocity  V  of  the  storm. 

Simplified,  steady  state,  one-dimensional,  analytical  models  with 
both  uniform  depth  and  uniform  slope  are  presented  with  solutions  in 
graphical  form.  Also  presented  is  an  analytical  treatment  of  a  one- 
dimensional  uniform  depth  shelf  subjected  to  a  triangular  shear  stress 
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distribution  moving  at  constant  velocity  C.  Dimensionless  maximum  re- 
sponse to  the  latter  treatment  is  given  as  a  function  of  dimensionless 
storm  v/idth  A*  and  velocity  c*.  These  solutions  allow  simplified  calcu- 
lations and  verification  of  the  numerical  methods. 

One-and  two-dimensional  models  are  developed  including  the  finite 
difference  equations  as  used  in  the  computer  simulation.  Several 
techniques  are  presented  for  verifying  the  results  of  the  numerical  cal- 
culations and  of  determining  the  effect  of  the  use  of  discrete  elements. 

Hurricane  Camille  of  August  17-22,1969  was  chosen  as  a  test  case 
for  the  numerical  techniques.  Two-dimensional  plots  of  surge  elevation 
at  various  times  are  included  as  well  as  the  surge  maximum  at  each  point 
along  the  coast;  the  latter  being  compared  to  high  water  marks  compiled 
by  the  U.S.  Army  Corps  of  Engineers.  The  sensitivity  of  the  surge  calcu- 
lations to  bottom  friction  coefficient  f,  and  to  starting  point  of  the 
storm  X  ,  was  investigated.  A  comparison  of  a  measured  tide  curve  to  the 
numerical  calculations  is  made  for  the  data  from  Dauphin  Island,  Alabama. 
The  currents  calculated  at  Eglin  AFB ,  Florida  are  compared  to  the  data 
taken  by  Murray  [16].  Techniques  to  be  used  for  field  measurements  of 
surges  and  storm  driven  currents  are  discussed. 

The  further  development  of  models  such  as  the  one  presented  will 
rely  on  a  better  understanding  of  the  physical  processes  involved. 
Interaction  between  the  wind  and  the  sea  surface  is  not  well  quantified. 
The  use  of  the  quadratic  bottom  friction,  predicated  on  steady  open  chan- 
nel flow,  can  probably  be  improved  upon.  The  availability  of  field 
measurements,  especially  during  hurricane  conditions  would  add  to  a  bet- 
ter understanding  of  the  subject. 


APPENDIX  A 
DERIVATION  OF  GOVERNING  EQUATIONS  FOR  SURGE  HYDRODYNAMICS 

In  this  appendix  the  equations  of  momentum  and  continuity  will  be 
developed  for  the  general  case  of  long  period  flows  in  estuaries,  bays, 
and  coastal  regions.  The  detailed  velocity  variations  over  depth  are  not 
considered  to  be  necessary  for  the  purpose  of  this  study.  The  problem 
is  therefore  formulated,  in  two  dimensions,  in  terms  of  the  vertically 
integrated  equations  of  motion  and  continuity  as  explained  in  Section  2-1 
of  the  main  text.  Referring  to  the  definition  sketch  Figure  2-1,  the 
equations  of  motion  are  Equations  (2-3)  of  the  main  text. 

t-lf-t-i-Msin,)v^-if.i(!^.!;i..!^^ 

lw^^|w   aw^^^  _1^_    l^!lk^1k^!lk)  (2.3c) 

at   "ax  ^dy      ^zz  p  3z   y   p  ^9x    8y    dz     '  ^^  ^^' 

The  vertical  accelerations  are  considered  to  be  small  as  are  stresses 
acting  on  any  vertical  plane. 

The  differential  form  of  the  equation  of  continuity  for  an  imcom- 
pressible  fluid  is: 

9x   ay   si  "  0  (2-3d) 

In  developing  the  vertically  integrated  equations  of  motion  it  will  be 
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necessary  to  utilize  Leibnitz'  rule: 


f(x,z,t)  dz 


h(x) 


g(x) 


|^(x,z,t)  dz  +  f(x,h(x),t)  1^ 


-f(x.g(x),t)  If 


Specifically  for  the  first  term  of  Equation  (2-3a)  we  have,  referring  to 
Figure  2-1  for  the  notation: 


"nCx.y.t) 


'-h(x,y) 


where  u_,  ||^ 


at 


dz  =  -^ 


^n(x,y,t) 


■h(x.y) 


^^^-^f^^^-hH^ 


0  since  h  is  not  a  function  of  time. 


To  render  Equations  (2-3a)  and  (2-3b)  into  better  forms  for  inte- 
gration, the  following  combinations  of  the  equations  of  motion  and  con- 
tinuity are  formed 


Eq.(2-3a)  +  u  •  [Eq.  (2-3d)] 
Eq.(2-3b)  +  V  •  [Eq.  (2-3d)] 


(A-1) 
(A-2) 


Rscalling  that  only  the  stresses  occurring  on  planes  perpendicular  to 
the  z  coordinate  will  be  considered  and  carrying  out  the  combinations 
indicated  by  Equations  (A-1)  and  (A-2)  the  resulting  equations  are: 


3U  .  3u2    3(uv)  .  9(uw)    r,    f     ■       \ 

at  ^37  +i-^r-^+^^-  2a.(sin*)v  = 


p  3X    p  9Z 


(A-3) 
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3t'. 


x^ 


p  ay   8z 


(A-4) 


Integrating  Equation  (A-3)  over  depth,  employing  Leibnitz'  rule,  and 
grouping  terms  leads  to  the  following: 


a_ 
3t 


r^ 


-h 


udz  + 


3X 


A 


-h 


A 


U2  d2  +  -^ 

3y 


/"n 


uv  dz 
■h         y 


2u)(sin4.)  V  dz 


-h 


n   3t    3x    ay    -"n 


+  u  ,  [  u^^^^+  v^-^^-  w  ] 
-h  ■-   ax      ay      -■ 


A 


1 

p 


J 


-h 


^     X      X 


(A-5) 


To  further  reduce  Equation  (^-S)   the  kinematic  free  surface  and  the 
bottom  boundary  conditions,  Equations  (A-6)  and  (A-7)  are  employed. 


an  .   an  ,   an 
at    ax    ay 


at  z 


(A-6) 


ax 


+  V 


a(-h) 
ay 


w 


at  z  =  -h 


(A-7) 


The  assumptions  stated  earlier  that  the  vertical  acceleration  is 
negligible  and  that  only  stress  on  horizontal  planes  is  of  importance, 
imply  that  the  motion  is  hydrostatic;  integrating  Equation  (2-3c)  over 
depth  yields  . 


p(x.y,z)  =  p^  (x,y)  +  Y(n(x,y)  -  z) 


(A-8) 


where  y  -■   pg  is  the  specific  weight  of  water  and  p  (x,y)  is  the  distri 
bution  of  (barometric)  pressure  on  the  free  surface. 
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We  now  have  the  necessary  equations  to  further  reduce  Equation 
(A-5).  Combining  Equations  (A-6),  (A-7)  and  (A-8)  with  Equation  (A-5) 
yields 

'n  f  n 

„2  dz  .  1^  I    uv  dz  =  -  1  (h.n)  ^   (A-9) 


u  dz  + 


at 


-h 


A        A 


The  form  for  the  equation  of  motion  in  the  y-di recti  on  is; 


9_ 
8t 


J 


^'^'k 


uv  dz  + 


3y 


-h 


J 


v2  dz  =  -  -  (h+n)  -P-      (A-10) 

p     3y 


■g(h+n)  ■^+r(T   -  T  ,  ) 


9y   P  ^'n^   '-h^ 


Integration  of  the  continuity  equation  over  depth  and  employing  Leibnitz' 
rule,  yields 

fn 


a_ 


u  dz  +  f- 
9y 


-h 


^^^-["lf*"t-"j„ 


-h 


.Cu|i^.v|i^-w].,  =0 


(A-ll) 


3X 


Combining  Equations  (A-6),  (A-7),  and  (A-ll)  leads  to  the  vertically 
integrated  equation  of  continuity.  Equation  (A-12). 
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In  +  i_ 
at   3x 


f 


u  dz  + 


3y 


V  dz  =  0 


(A-12) 


-h 


The  governing  equations  will  be  more  convenient  for  later  use  if 
the  following  definitions  are  made  and  the  form  of  the  equations  modi- 
fied accordingly: 


u  dz 


r  n 


y., 


V  dz 


D(x,y,t)  =   h(x,y)  +  n(x,y,t) 

The  vertically  integrated  equations  of  motion  and  continuity  are  sum- 
marized as: 


3q 


X  ^  9 


9t     3X 


u2  dz  +  1^ 


D  ^P 
uv  dz  -  2a)(sin(f>)qy  "  "  ^  9x^  "  ^D 


9n 
3x 


+  -  (t   -  t  u  ) 
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at   ax 


uv  dz  + 


ay 


V''-  dz  +  2w(sin(}))q^ 


p  ay   ^  ay 


*  F  K  ■  '-h  > 


!^.!il,ln 


ax 


9y 


at 


To  express  these  equations  completely  in  finite-difference  form, 
further  approximations  must  be  made  as  discussed  in  Appendix  B. 


APPENDIX  B 

TRANSFORMATION  OF  GOVERNING  DIFFERENTIAL  EQUATIONS 
INTO  FINITE  DIFFERENCE  FORM 

1.  Transformation  Without  Convective  Terms 


Transforming  the  vertically  integrated  equation  of  motion  into  a 
finite  difference  form  is  straightforward  if  the  nonlinear  friction  and 
convective  terms  are  neglected.  If  they  are  not  to  be  neglected  then 
some  form  of  approximations  must  be  made  in  addition  to  the  approxima- 
tions introduced  by  discretization.  Writing  the  vertically  integrated 
equations  of  motion  and  continuity  and  including  the  convective  terms 
yields  Equations  (B-la),  (B-lb),  and  (B-lc). 


at   3x 


dz  +  — 
9y 


uv  dz  -  2  CL)(sinA)q  +  gD  |^ 


-h 


-D  ^n^l      ^1^'^) 
p  3x    p  ""n^     8D^ 


(B-la) 


3q 


X+  3_ 


3t    3X 


J 


-h 


r 


uv  dz  .  I3; 


J 


9n  _ 


v2  dz  +  2  a3(sin4.)q^  '^  ^^  ^ 


p  3y   p  ■'n     80^ 


(B-lb) 
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in  +   —JL  +   _X  :=    0 
3t        9X  3y 

''  n  fn. 


(B-lc) 


where,  as  before,  q 


J 


u  dz,  q  =1   V  dz,  and  where  the 


Darcy-Weisbach  friction  coefficient,  f,  is  used,  to  represent  the 
bottom  friction. 

The  numerical  computations  presented  herein  v/ere  carried  out 
without  the  convective  terms  represented  by  the  integrals  in  Equations 
(B-la)  and  (B-lb),  These  equations  then  assume  the  form  of  Equations 
(3-6)  in  the  text. 


9t 


-  2co(sin({.)qy  +  gD  ^^ 


9n  _  -D 


9p 


p  3X 


n  +  i 


8D^ 


{3-6a) 


3q 


^+  2a)(sin(^)q^  +  gD  |^=  - 


9y 


p  ay    p  ^y^   ~80^ 


(3-6b) 


apv   S^v 


9t   3X 


9y 


(3-6c) 


To  transform  Equation  (3-6a)  into  finite  difference  form  integrate 
from  t  to  t  +  At.  Using  the  indicia!  notation  described  in  Chapter 
3  of  the  text,  the  result  of  that  integration  is  approximately 
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q^^Ci  J.n+1)  -  q^(i,j,n)  =  2  uCsiricj,)   [ 


n+l 


qy(i,J,C)   dc] 


+  ?7t[  D(i.j,n+l/2)  +  D(i-l,j.n+l/2)]    [n(i-l  .J,n+l/2)   -  n(i  ,j»n+l/2)] 


2ax 


At 


+  f^[  D(i,j.n+l/2)  +  D(i-l,j.n+l/2)]   [p^(i-l  .j  ,n+l/2)   -  p^(i  ,j  ,n+l/2)] 


+  ^l\     (i.J.n+1/2)  +  T^     (i-l.j.n+1/2)] 
n  n 


n+l 


f/2 

LDC1,J,n+l/2)  +  D(i-1  J,n+l/2)]2 


|q(iJ,C)I   qJiJ,c)  d?   (B-2) 


Where  i  is  a  dumiriy  variable  and  q   (i,j,c)   is  defined  in   Figure  3-4. 

The  value  of  D  at  the  location  of  q   (i,j,n)   is  assumed  to  be  the  average 

of  D(i,j,n)  and  D(i-l,j,n)  and  furthermore  the  value  at  time   (n+1/2)   is 

assumed  to  be  adequately  represented  by  the  average  value  from  time   (n) 

to  time   (n+l).     To  evaluate  the  integral   terms  in  Equation  B-2     the 

assumption  is  first  made  that  the  parameters  q   ,  q     are  slowly  varying 

X   y 

with  respect  to  the  time  scale  At. 


At 


t  +  At 


qy(i.j,c)  dc 


is  then  approximated  by  q  (i,j,n)At,  remembering  that  q  (i  ,j,n+l/2)At 
would  be  a  better  approximation  but  that  the  value  at  (n+1/2)  is  not 
available.  An  iterative  technique  could  be  used  to  obtain  q'y(i  J  ,n+l/2) 
but  the  magnitude  of  the  error  is  small  since  F  is  seen  to  change  slow- 
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ly  with  respect  to  At  (see  Figure  4-10),  and  cost  in  increased  computer 

time  is  large  so  that  a  change  to  an  implicit  scheme  is  not  warranted 
at  the  present  time. 

A  more  critical  problem  is  that  of  placing  the 


n+1 


iq(ij,c)|  q^Ci.J.S)  dc 


term  into  finite -difference  form.     To  maintain  an  explicit  form  the 
following  approximation  is  made: 


n+1 


|q(i,j,c)|   qxd'.J.S)   dc  =    |q(i,j,n)|   q   (i,j,n+l)  At 


To  show  that  this  is  correct  to  second  order,  it  is  sufficient  to 
demonstrate  that  to  second  order  in  At 


/- 


t^+At 


q2dt  =  qCt^)  q(tQ+At) 


Expanding  about  q(t  ) 


to*" 


J  t. 


ft  *  1^  t'-'o'  *  0  C-'o'']'  ■" 


ft„+At 


[,2.2q  |a  (t-tj  *   ,  0  (t.tj^ttia-j^t-t^jViJdt 
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Making  the  change  of  variable  t'  =  t  -  t 


At 


h^-^^'ift'^^^t'^.^ClfjS'^ldf 


=  q  [q  +  2||-  At  +  0(At2)]At 


Now  changing  back  to  t  =  t'  +  t  and  expressing  the  result  in  indicial 
notation  yields  to  second  order  in  At, 

=  qd.J.n)  [q(i,j,n+l)]At 

The  desired  result  with  the  introduction  of  the  x-component  and 
absolute  value  is 

|qCi.J.n)I   q^(i,j,n+l)At 

Introducing  this  into  Equation  (B-2a)  and  solving  for  q  (i,j,n+l) 
yields  Equation  CB-3a) 

q^Ci.J.n+l)  =  Y~   C^x^^'^'"^  ^  2u(sin({))q  (i  ,j,n)At 

A 

+  ^   {D(i,j,n+l/2)+D(i-l,j,n+l/2)}{^P^(i-l,j,n+l/2)-p^(i,j,n+l/2)] 

+  g  [n(i-l,j,n+l/2)  -  n(i  ,j  ,n+l/2)]} 

+  1^  {t   (i,j,n+l/2)  +  t   (i-l,j,n+l/2)}]  (B-3a) 

^x  ^x 
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with 

F^  =  1  +^{[q^(iJ.n)]2  +  [q^(i  .j,n)]2}  .  (B-3b) 

[D(i-l,j,n+l/2)  +  D(i,j,n+l/2)]"^ 

These  are  Equations  (3-7a)  and  (3-7b)  of  the  main  text  respectively. 
It  can  be  seen  that  forward  differences  in  time  have  been  used 
throughout.  Application  of  this  procedure  to  the  y-equation  and  the 
continuity  equation  is  identical  in  concept  to  the  above  and  the 
equations  in  finite  difference  form  are  given  by  Equations  (3-8a), 
(3-8b),  and  (3-9)  of  the  main  text. 
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2.     Effect  of  Convective  Terms 
The  convective  terms  for  the  x-equation  are:     |—  /         u^  dz  and 

oX    J 

-h 


J- I       uv  dz.     Since  the  form  of  the  velocity  profile  is  an  unknov;n, 
-h 

evaluation  of  the  convective  terms  is  difficult.     An  approximation 

could  be  made  by  assuming  that  the  flow  is  highly  turbulent   (   a 

physically  reasonable  condition),  and  similar  in  form  to  that  of  open 

channel   flow,   that  is  steady  state.     This  would  then  provide  a  velocity 

distribution  that  is  logarithmic  in  form  and  reasonably  uniform  above 

a  distance  0.2  x  depth.     Making  the  more  restrictive  assumption  that  the 

velocity  distribution  in  the  vertical   is  uniform  enables  one  to  evaluate 

the  integral   terms. 


^  u2  dz  = 


8X    L    ^^    -I         9x    I-    D    -^ 


uv  dz 


^  [  Duv  ] 


9_r  iii 
ay  •-    D     -■ 


where  q^  =  Du  and  q     =  Dv  as  a  result  of  the  uniform  velocity  distribu- 
tion. 

If  one  continues  at  this  point  as   in  part  1   of  this  appendix  and 


assumes 
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q^(n+l/2)  q^(n+l/2)  =  q^(n+l)  q^Cn) 

an  implicit  scheme  would  still  be  necessary  because  of  the  required 

spatial  derivatives  in  x  and  y.  Keeping  in  mind  the  magnitude  of  the 

error  in  part  1  of  this  appendix  it  is  reasonable  to  approximate  q  ^, 

q  q  ,  and  q  ^  by  their  values  at  one  time  step  earlier,  this  is  es- 
X  y      y 

pecially  true  near  the  time  of  peak  surge  and  currents  where 


9qv      9q. 

"ion 


^   0  and  TT^ -»■  0.  The  x-equation  of  motion  then  becomes  Equati( 


at       at 
(B-4) 

qx(i'J'"+'')  "  -f-   [Px^^'J'"^)  ^   2co(sin(j.)q^(i,j,n)  At 

+  H^  {D(i-l,j,n+l/2)  +  D(i,j,n+l/2)}  {  ^  [p^Ci-l  ,j  ,n+l/2)-p^(i  ,j  ,n+l/2)] 
+  g  [  n(i-l,j,n+l/2)  -  n(i,J,n+l/2)  ]} 

+  fj  {t  (i-l,j,n+l/2)  +  T   (i,j,n+l/2)} 

^   ^x  X 

-  ^  [  D(i-l.j,n+l/2)  +  D(i,j.n+l/2)]'^  [q/(i+l,j.n)  -  q/(i-l,j,n)] 

-  f^   [  D(i-l.j,n+l/2)  +  D(i,j,n+l/2)]'^[q^(i,j+l.n)  ^^(i.j+l.n) 

-qx(i.J-''.i)  qy(i.a-"i>n)]]  (b-4) 

where,  as  before 

F^  =  1+  ^  {  [qx(i.j>n)]2  +  [?y(i,j,n)]2  }l/2  . 

[  D(i-l,j,n+l/2)  +  D(i,j,n+l/2)]"^ 

The  convective  terms  for  the  y-equation  of  motion  are 
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uv  dz  and 


3y 


v^dz. 


J 


applying  the  approximations  used  for  the  x-equation,  the  y-equation  can 
be  written  as  Equation  (B-5). 

qy(ij,n+l)  =  p-  [  q^(i,j,n)   -  2a)  sin<j)  q'^(i,j,n)At 


At 


+  2^3^  {D(i,j,n+l/2)  +  D(i,j-l.n+l/2)} 


{  ^[p^   (i,J-l,n+l/2)   -  p^(i,j,n+l/2)   ] 


+  g  [n(i,J-l,n+l/2)   -  n(i  ,j,n+l/2)   ]   } 


+  0  {t^   (i,j,n+l/2)  +  T^   (i,j-l,n+l/2)   } 


At 


4^  [D(i,j,n+l/2)  +  D(i,j-l,n+l/2)]   [q   (i+l,j,n)  q^(i+l,j,n) 


-qy(i-l.j,n)   q^(i-l,j,n)] 


At 


4^y  [D(i,j,n+l/2)  +  D(i.j-l,n+l/2)]   [q^Mi  .J+1  ,n)   -  q^2(i  ,j_i  ^p)]] 


(B-5) 


where 


fAt 


Fy  =   1   +^{[q^(iJ,n)]2  +  q   (i.j.n)]^} 


-1/2 


[D(i,j,n+l/2)  +  D(i,j-l,n+l/2)]" 

A  better  smoothing  of  D  as   in  Section  3-b  of  the  text  may  be  desirable 
but  the  simplicity  of  the  averages  used  in  Equations   (B-4)  and   (B-5) 
dictate  that  they  be  used  since  the  goal   here  is  to  obtain  a  numerical 
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estimate  of  the  changes  due  to  the  inclusion  of  the  convective  terms. 

Introduction  of  these  forms  for  the  convective  terms  into  the 
numerical  model  changed  the  calculated  maximum  surge  by  about  1.2% 
and  changes  the  surge  height  at  any  grid  element  adjacent  to  the  coast 
by  less  than  5%.  Because  of  the  complexity  and  their  relatively  small 
effect,  the  convective  terms  were  not  further  considered  in  this  study. 
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